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Abstract 


Purpose  of  the  present  investigation  was  to  assess  the  feasibility  of  simulating  and  study¬ 
ing  coherent  structures  in  turbulent  shear  layers,  making  use  of  Large  Eddy  Simulations 
(LES). 

The  current  investigation  has  been  performed  across  a  rather  wide  range  of  simulated 
flows  and  LES  implementations,  starting  from  an  established  and  validated  structured  LES 
code  based  on  finite  differences  and  cartesian  grids  ending  to  an  unstructured  LES  code  un¬ 
der  development.  The  use  of  the  unstructured  approach  was  considered  necessary  whenever 
complex  geometries  have  to  be  analyzed,  which  cannot  be  handled  by  the  structured  code, 
regardless  its  advanced  features  (like  multi-domain  approach).  The  unstructured  LES  code 
was  developed  in  the  frame  of  current  GRANT,  and  it  has  undergone  validation.  The  struc¬ 
tured  LES  code,  which  had  already  validated  over  a  wide  range  of  flow  configurations,  was 
used  to  analyze  the  wake  behind  an  obstacle,  attempt  which  was  not  made  before.  Extensive 
studies  were  needed  in  the  attempt  to  validate  to  code  for  this  case. 

Given  the  wide  scope  of  present  report  and  the  mass  of  data  available  from  different 
LES  codes,  flows  and  structures,  the  work  has  been  spit  between  two  work  groups,  the  one 
devoted  to  the  structured  LES  and  the  one  devoted  to  the  unstructured  LES.  In  this  frame 
it  was  found  more  practical  to  sub-divide  the  full  report  in  three  volumes: 

Volume  I  Investigation  of  Coherent  Structures  in  an  Attached  Shear  Layer 
Volume  II  Wake  around  square  Cylinder 
Volume  III  Unstructured  LES 

The  different  volumes  are  not  completely  independent:  Volume  I,  for  example,  in  the 
description  of  the  structured  LES  code,  introduces  concepts  and  algorithms  which  will  be 
applied  in  computations  that  will  be  presented  only  in  Volume  II. 

The  subdivision  of  present  report  in  different  volumes  was,  indeed,  adopted  to  simplify 
the  organize  and  present  the  collected  data,  and  also  to  take  into  account  the  fact  that 
different  groups  worked  on  different  aspects  of  the  grant,  and  their  respective  work  should 
be  easily  recognizable. 

Volume  I  of  present  report  has  described  the  general  background,  the  main  development 
and  the  results  obtained  for  wall  bounded  turbulent  flow,  particularly,  for  the  case  of  turbulent 
plane  channel  flow  at  equilibrium. 

Volume  II  is  devoted  to  the  study  of  the  flow  around  an  obstacle,  specifically,  a  cylinder 
with  square  cross  section  and  infinite  length  in  span- wise  direction.  This  test  case  belongs 
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to  a  completely  different  class  of  flows,  respect  to  the  one  analyzed  in  Volume  I,  and  it  is  a 
severe  test  case  for  the  LES  code. 

Most  of  the  advanced  features  of  the  code  were  exploited  to  achieve  this  task,  different 
simulations  were  performed,  in  order  to  collect,  whenever  possible,  sufficient  data  to  draw 
conclusions  and  remarks. 

Differently  from  the  case  of  channel  flow,  where  the  code  has  been  extensively  validated 
and  analyzed  in  detail,  down  to  the  level  of  turbulent  kinetic  energy  budget,  the  present  flow 
constitutes  a  still  active  area  of  basic  research. 

The  higher  complexity  of  the  flow,  the  absence  of  a  second  direction  of  homogeneity, 
beside  the  span-wise  direction,  the  existence  of  a  dominant  frequency  in  the  flow,  the  necessity 
of  employing  upwinding  discretization  for  advection  terms  and  other  numerical  difficulties 
have  led  to  such  complexities  of  execution  that  only  the  very  basic  set  of  statistics  (statistical 
moments  of  first  and  second  order)  were  obtained.  At  the  same  time,  limits  and  draw-backs 
for  current  LES  code  implementation  have  been  outlined. 

Most  of  the  efforts  were  devoted  to  reach  a  sufficient  similarity  in  behavior  between  LES 
simulation  and  reference  experimental  data;  this  attempt  has  led  to  build  grids  not  perfectly 
suited  to  analyze  structures  applying  the  criteria  introduced  in  Volume  I,  which  are  based 
on  operators  working  on  the  tensor  of  gradient  of  velocity. 

This  fact,  together  with  intrinsic  difficulty  of  dealing  with  the  presence  of  a  wide  range 
of  structures  of  different  sizes  and  shapes,  has  limited  the  possibility  to  obtain  a  through  and 
complete  description  of  vortical  behavior  for  this  class  of  flow. 

However,  it  was  possible  to  find  useful  global  statistics  for  organized  turbulence  and  to 
show  that  a  very  important  part  of  the  entire  turbulent  field  lies  within  very  active  structures, 
which  cover  a  very  , small  part  of  the  entire  flow  field.  It  can  be  moreover  suggested  that  these 
characteristics  of  these  structures  are  the  ones  which  define  and  control  the  wake  behavior, 
and  their  identification  could  be  key  of  understand  and  possibly  identify  and  recognize  specific 

wakes. 

Possible  avenues  for  future  investigations  starting  from  present  study  are  also  suggested. 
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1.  INTRODUCTION 


Since  it  was  first  recognized,  turbulence  has  defied  any  complete  definition: 

Observe  the  motion  of  the  water  surface,  which  resembles  that  of  hair,  that  has  two 
motions:  one  due  to  the  weight  of  the  shaft,  the  other  to  the  shape  of  the  curls;  thus,  water 
has  eddying  motions,  one  part  of  which  is  due  to  the  principal  current,  the  other  to  the  random 
and  reverse  motion.  (Leonardo  da  Vinci  1510) 

Under  the  practical  point  of  view,  the  most  visible  and  striking  feature  of  turbulent  flows 
would  appear  their  exhibited  unsteadiness  and  irregularity,  leading  to  consider  the  details  of 
their  behavior  random  and  unpredictable. 

However,  starting  from  the  experiments  of  Brown  and  Roshko  in  1975  (Brown  and 
Roshko,  1975),  the  idea  that  there  could  be  an  underlying  order  in  the  chaos  of  the  turbu¬ 
lence  has  gained  momentum;  the  difficulty  in  the  matter  lying  with  its  mix  of  deterministic 
organization  together  with  chaotic  features. 

...frustration  results  largely  from  the  mixture  of  chaos  and  order  and  the  wide  range  of 
length  and  time  scales  that  turbulent  flows  possess.  (Mom  and  Mahesh,  1998) 

The  concept  of  an  underlying  deterministic  structure  of  turbulent  flow  has  been  expressed 
by  the  definition  of  coherent  structure,  originally  given  by  Hussain  (Hussain,  1983,  1986). 


A  coherent  structure  is  a  connected  turbulent  fluid  mass  with  instantaneously 
phase  correlated  vorticity  over  its  spatial  extend.  ...We  choose  to  designate  the 
instantaneously  space  and  phase- correlated  vorticity  as  coherent  vorhcvtij.  Thus, 
coherent  vorticity  is  the  primary  identifier  of  coherent  structures,  which  have  dis¬ 
tinct  boundaries  and  independent  territories. 


Therefore,  the  understanding  of  the  physics  of  the  coherent  structures  and  their  inter¬ 
action  with  the  background  motion  is  of  paramount  interest  in  turbulence  research.  Indeed, 
recently  had  been  suggested  by  Farge  (Farge  et  al.,  1996)  that  separating  and  understanding 
the  different  physical  behavior  of  coherent  structures  and  un-coherent  background  might  be 
the  key  to  understand  and  predict  the  turbulence  behavior. 

Most  previous  studies  on  the  topic  followed  the  experimental  approach:  detection  of 
coherent  structures  and  deduction  of  their  behavior  alongside  their  physics  are  obtained 
by  velocity  measurements  of  the  flow,  applying  point  measurements  as  such  Laser  Doppler 
Velocimetry  (LDV)  (Lyn  et  ah,  1995);  multi-point  point  measurements  such  as  a  rack  of 
Hot  Wire  (HW)  probes  (Hussain  and  Hayakawa,  1987);  in  plane  measurements  as  such  as 
Particle  Image  Velocimetry  (PIV)  (Scarano  et  ah,  1999).  Of  course,  in  all  the  cases  previously 
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considered,  it  is  impossible  to  have  simultaneous  three-dimensional  data  to  take  into  account 
the  inherited  three-dimensional  nature  of  turbulent  processes. 

The  fundamental  advantage  of  the  numerical  approach,  in  the  present  case,  is  its  capa¬ 
bility  to  produce  instantaneous,  three-dimensional  flow-fields.  The  tremendous  potential  of 
this  approach  was  demonstrate  by  pioneering  work  of  Robinson  k  al.  (Robinson  et  al.,  1990) 
who  visualized  the  coherent  structures  in  boundary  layer  flow  over  a  flat  plate. 

This  study  was  performed  using  flow-fields  created  by  Direct  Numerical  Simulation 
(DNS)  (Spalart,  1988)  where  all  the  turbulent  scales  down  to  the  Kolmogorov  scale  are 
resolved. 

While  the  advantages  of  DNS  are  evident  and  well  know  (Moin  and  Mahesh,  1998),  its 
numerical  costs  are  likely  to  confine  its  applications  to  Reynolds  numbers  (Re)  and  simple 
geometry  for  next  future. 

Simulation  of  flows  at  higher  Re  and  complex  geometry  could  be  feasible  applying  Large 
Eddy  Simulation  (LES)  (Rodi  et  al.,  1997;  Piomelli,  1997). 

The  basis  of  LES  approach  is  the  application  of  a  filter  operation  to  the  Navier-Stokes 
equations  to  obtain  the  corresponding  transport  equation  for  filtered  scales:  scales  larger  than 
width  of  filter  (grid  scales)  are  resolved  by  numerical  solution.  Scales  than  smaller  width  of 
filter  (sub-grid  scales)  are  cut-off  by  filtering  operation;  effect  sub-grid  scales  on  the  grid  ones 
takes  the  form  of  stresses  (sub  grid  stresses:  SGS)  and  are  modeled  with  an  ad  hoc  model 
(sub-grid  model).  Under  the  physical  point  of  view,  the  concept  is  to  reproduce  the  different 
behavior  of  large  eddies  and  small  eddies  in  the  transfer  of  turbulent  kinetic  energy  (Tennekes 
and  Lumley,  1983):  assuming  the  spectrum  of  energy  to  have  an  equilibrium  range,  it  can  be 
supposed  that  smaller  scales  in  this  range  are  nearly  isotropic  in  their  behavior,  and  therefore 
easier  to  be  modeled  than  the  full  turbulent  field. 

Available  understanding  of  coherent  structures  shows  that  their  size  is  considerably  larger 
than  the  turbulent  macro-scale,  making  permissible  the  application  of  LES  to  their  simulation 
and  study,  with  the  advantage  to  be  able  to  study  their  behaviour  at  higher  Re  and  complex 
geometry. 

Purpose  of  the  present  study  is  to  validate  the  use  of  LES  to  the  study  of  coherent 
structures  in  shear  layer  and  apply  it  to  different  flow  cases  at  high  Re. 

As  previous  studies  have  shown  that  structures  characteristics  are  quite  different  for 
different  classes  of  flow,  notably  between  wall-bounded  and  free  flows,  three  different  cases 
will  considered,  namely: 

1.  Attached  wall-flows,  for  which  plane  channel  at  equilibrium  will  be  adopted  as  test  case 

2.  Free  flows,  for  which  case  plane  mixing  layer  will  be  studied 
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3.  Separated  flows;  which  will  be  represented  by  flow  around  a  bluff  body 

All  these  three  types  of  flow  are  investigated  in  order  to  find  out  the  general  and  universal 
characteristics  of  the  coherent  structures. 

The  initial  point  of  the  present  work  is  the  Large  Eddy  Simulation  code  developed  at 
the  von  Karman  Institute  (Simons,  2000). 

Its  fundamental  features  are: 

•  Finite  differences 

•  Cartesian  coordinates 

•  Staggered  Grid  Approach 

•  Second  order  accurate  in  space  and  time 

•  Filtering  implicitly  performed  from  finite  differences 

•  Multi-Domain  approach  for  not  trivial  geometries 

In  order  to  apply  this  code  the  study  of  the  coherent  structures  the  main  task  to  be 
performed  will  the  development  of: 

•  Algorithms  for  detection  for  the  coherent  structures  in  three-dimensional  fields 

•  Statistical  tool  for  analysis  of  the  structures 

These  capabilities  are  implemented  in  a  code  treating  structures  defined  on  structured, 
Cartesian  grids. 

In  order  to  assess  the  capability  of  LES  to  the  study  of  coherent  structures  in  general 
flows  which  cannot  well  be  represented  over  Cartesian  grids,  a  code  which  solves  laxge-eddy 
equations  in  un-structured  frame  has  been  developed. 

Present  report  is  the  subdivided  in  3  volumes,  each  for  the  basic  topics  of  the  study. 

Volume  I  Is  devoted  to  the  discussion  of  the  structured  code,  the  development  of  algorithm  for 
detection  and  classification  of  the  coherent  structures  and  assess  the  result  of  its  appli¬ 
cation  to  simulation  of  attached  wall  flow. 

Volume  II  Covers  the  simulation  of  flow  and  coherent  structures  around  the  cylinder  of  square 
cross-section,  and  assess  the  performances  of  the  present  algorithms  for  wake  flows. 
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Volume  III  Discusses  the  development  and  the  performances  of  unstructured  LES  code. 

Present  Volume  I,  then,  discusses  basic  elements  of  phomenology  of  coherent  structures 
and  all  mathematical  and  statistical  tools  developed  for  they  analysis.  Performances  of  these 
algorithms,  sensitivity  of  simulation  to  different  physical  and  computational  parameters  and 
their  influence  on  the  definition  the  coherent  structures  are  assessed  in  the  case  of  plane 
turbulent  channel  at  equilibrium. 

The  Volume  is  organized  as: 

•  Theory  of  LES  is  briefly  discussed  in  §  2 

•  Structured  LES  solver  is  presented  in  §  3 

•  Development  of  code  for  detection  and  analysis  of  coherent  structures  is  presented  in 
§4 

•  Study  of  coherent  structures  in  attached  wall-flow  is  presented  in  §  5  and  §  6 

•  Achieved  conclusions  are  presented  in  §  7 
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2.  LARGE-EDDY  SIMULATION  CONCEPT 


2.1  Governing  Equations 

The  Navier-Stokes  equations  provide  a  simplified,  coarse-grained  description  of  a  Newto¬ 
nian  fluid  flow  and,  despite  the  enormous  simplification  they  represent  respect  of  the  molec¬ 
ular  dynamics,  contain  all  the  relevant  physics  of  the  fluid  at  scales  much  larger  than  the 
mean  free  path  of  the  molecules.  As,  even  the  smaller  scales  occurring  in  a  turbulent  flow 
are  ordinarily  far  larger  than  any  molecular  length  scale  (Tennekes  and  Lumley,  1983)  they 
allow  to  study  and  simulate  and  all  the  features  of  turbulent  processes. 

The  present  study  is  restricted  to  flows  at  low  Mach  numbers,  where  density  and  molecu¬ 
lar  viscosity  can  be  assumed  constant.  In  this  case,  the  governing  equations  reduce  to  the  well 
know  incompressible  Navier-Stokes  equations  for  the  conservation  of  mass  and  momentum, 
whose  adimensional  formulation  is: 


dui 

dxi 

dui  duiUj 
dt  +  dxj 


0 


dxi  Re  dx{dxj 


(2.1) 

(2.2) 


where  length:  X{ ,  velocity:  U{  and  pressure:  P  are  adimensionalized  against  a  reference  length 
Sref  and  velocity  uref  (superscript  (*)  denotes  dimensional  quantities): 


Ui  — 


-  r-  )  **1 

STef  uTef 


U*  t  Ursf 

Sr 


P  = 


?ref  P^ref 

and  Re  is  the  Reynolds  number  based  on  the  above  reference  quantities. 


(2.3) 


Re  = 


firefflref 


(2.4) 


2.2  Principles  of  Direct  and  Large-Eddy  Simulation 

Turbulence  consists  of  a  continuous  spectrum  of  scales  ranging  in  size  from  large  scales  of 
the  same  order  of  magnitude  as  the  geometry  of  the  flow  to  smallest  scales  whose  correspond¬ 
ing  Re  is  of  the  order  of  1,  which  are  also  known  as  Kolmogorov  scale  (Tennekes  and  Lumley, 
1983;  Ferziger,  1997).  Kolmogorov  introduced  the  idea  of  an  energy  cascade  (Fig.  2.1)  which 
assumes  that  the  turbulent  energy  produced  by  the  largest  turbulent  structures  is  transferred, 
via  the  intermediate  eddies,  toward  the  smallest  eddies  where  this  it  is  dissipated  by  viscous 

effects. 
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Fig.  2.1:  Typical  spectrum  of  kinetic  energy  of  turbulent  flows 

Spectral  analysis  of  kinetic  energy  associated  with  turbulent  field  E,  versus  wave- number 
k  shows  that,  for  flows  characterized  by  developed  turbulence  ,  the  length  scales  can  be 
categorized  into  three  distinct  classes  (Tennekes  and  Lumley,  1983;  Ferziger,  1997): 


1.  “Large”  scales :  the  integral  length  scales  L  in  the  flow  are  non-universal,  highly 
anisotropic,  due  to  the  effect  of  the  geometry  on  the  flow.  The  wave  numbers  k  are  of 
the  order  of  27r/L. 

2.  “Medium”  scales:  scales  for  which  the  energy  spectrum  E(k)  becomes  more  universal, 
moderately  isotropic,  and  can  be  described  by  equation  2.5. 

E(k)=CK£2/3k-5^  (2.5) 

where  where  Ck  is  the  Kolmogorov  constant,  £  is  the  dissipation  rate. 

3.  “Small”  scales:  scales  for  which  the  energy  spectrum  E{k)  is  universal,  highly  isotropic, 
and  at  which  the  energy  is  finally  dissipated.  The  wave  number  is  of  the  order  k  = 
Re^L-K 


Simulation  of  turbulence  by  numerical  means  then  requires  to  solve  the  Navier-Stokes 
equations  marching  in  time  with  a  space  discretization  refined  enough  to  resolve  the  smallest 
scales.  Consequently  Direct  Numerical  Simulation  (DNS)  is  obviously  the  preferred  tool  for 
investigation  of  the  physics  of  turbulence  (Moin  and  Mahesh,  1998),  and,  in  present  case,  for 
the  study  of  the  phenomenology  of  coherent  structures. 

However,  its  application  is  and  is  likely  to  remain,  for  some  time,  confined  to  flow  at 
low  SR  and  simple  geometry,  in  view  of  the  fact  that  the  spatial  resolution  necessary  for  the 
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simulation  of  the  smallest  scales  increases  with  Re 9/4  while  the  corresponding  decreasing  on 
the  time  scales  leads  the  computer  time  of  simulations  to  increase  with  a  factor  of  Re  . 

Simulations  at  higher  Re  (typically  Re  >  104,  at  present  time)  could  be  reasonably 
feasible  only  applying  the  Large  Eddy  Simulation  (LES). 

This  modeling  concept  takes  advantage  of  the  different  behavior  of  turbulence  at  largest 
and  smallest  scales  (see  above)  to  separate  (Fig.  2.2)  the  large,  highly  anisotropic,  energy 
carrying  scales  which  are  directly  resolved  from  the  small,  nearly  isotropic,  energy  dissipating 
scales  which  are  modeled. 


grid-  scale  turbulence  turbulence 


Log(K) 


Fig.  2.2:  Separation  of  grid-scale  and  sub-grid  scale  turbulence 


Separation  of  the  different  scales  is  achieved  applying  a  filtering  operator,  which  takes 
the  form  a  space  filter  expressed  as  a  convolution  integral  (Germano,  1986,  1992). 


u,  (1, t)=  jj jf  <?  (f  -  a)  Ui  (6 1) d3i  (2.6) 


In  case  of  Finite  Difference  (FD)  and  Finite  Volume  discretization  the  well  known  ^top 
hat”  filter  (Piomelli,  1997) 


x-  f 


<  A/2 


otherwise 


(2-7) 


is  implicitly  applied  by  the  grid  selected  for  the  simulation,  size  of  filter  A  corresponding  to 
size  of  the  grid  A Xi  . 

Concerning  the  formulation  of  spatial  derivatives  of  filtered  quantities  it  must  remarked 
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that,  on  non-uniform  grids,  filtering  and  differentiation  do  commute  only  at  a  2nd  order 
commutative  error  (Ghosal  and  Moin,  1995),  which  should  be  taken  in  account,  in  case  of 
application  of  discretization  of  order  higher  than  2. 

Application  of  filtering  operation  to  the  incompressible  Navier-Stokes  equations,  yields 
the  LES  equations: 


dxi 

duj  dufilj  _  dP  1  d2Uj  |  - 
dt  +  dxj  ~  dxi  Redxjdxj  1 

where  the  over-bar  denotes  filtered  (resolved)  quantities.  As  resolved  and  un-resolved  scales 
are  defined  by  grid  size,  it  is  current  to  use  the  definition  of  “grid”  and  “sub-grid”  scales  for 
the  two  classes. 

To  express  and  solve  the  above  equation  in  term  of  U{  it  is  necessary  to  replace  the  term 
ufu~  by  the  one  UiUf,  the  difference  between  the  two  terms  takes  the  shape  of  a  tensor 

Tij  =  UiUj  —  UiUj  (2.10) 

which,  by  analogy  with  Reynolds  tensor,  is  named  “sub-grid  stress”  (SGS)  tensor. 

Equation  2.8  becomes 


(2.8) 

(2.9) 


dui  .  duiUi  dP  ,  1  d2Ui  dry  ,  T 

aT+  fef  =  ~d^P'Re8^r^fi 

since  each  tensor  can  be  decomposed  in  its  isotropic  and  anisotropic  components 


and  the  isotropic  part  of  the  SGS  tensor  can  be  added  to  the  pressure  obtaining  a  modified 
pressure  term 


Finally  2.8  becomes 


P  =  P  +  \rkk 


(2.13) 
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cm,. 

dt 


dui 

dxi 

=  0 

dr{a)  - 

-  *  +fi 

(2.14) 

duiUj 

dp  1  d2Ui 

(2.15) 

dxj 

~~  dxi  Re  dxjdxj 

axj 

(2.16) 

obviously,  r-f  must  be  modeled  (SGS  model)  to  close  the  problem. 

The  major  advantage  of  LES  concept  with  respect  to  Reynolds-averaged  Navier-Stokes 
(RANS)  approach  lies  on  the  fact  that,  while  the  latter  attempts  to  model  all  turbulent 
scales,  including  the  highly  anisotropic  scales  of  small  wave-numbers,  LES  models  only  the 
effect,  mostly  dissipative,  of  smaller  scales,  which  could  be  thought  isotropic,  at  least  m  firs 
approximation,  and  could,  therefore,  be  expected  to  be  easier  to  model. 

It  has  to  remarked  that,  refining  the  grid  (reducing  filter  size)  size  of  rjf  term  reduces 
itself  till  all  the  relevant  scales  are  captured  and  the  DNS  solution  is  recovered: 


Ui  =  Ui 


(2.17) 


2.3  SGS  Modeling 

As  stated  above,  the  r{j  term,  which  must  be  somehow  modeled,  represents  the  interac¬ 
tion  between  the  unresolved  small  scales  and  the  resolved  large  scales.  Different  models  are 
proposed  by  the  current  literature  and  is  reviewed,  among  others,  by  Piomelli  (1997),  for  t  e 
present  study,  the  attention  SGS  model  is  justified  by  the  need  to  verify  to  ensure  that  the 
simulated  flow-field  well  reproduces  the  statistics  of  the  turbulence  and  the  need  to  assess 
the  possible  influence  of  SGS  model  on  the  simulation  of  structures  themselves;  therefore 
following  discussion  is  confined  with  models  actually  applied  in  the  frame. 

SGS  models  can  be  sub-divided  in  two  classes  :  fully  dissipative  and  self-similar  models. 
2.3.1  Fully  dissipative  SGS  models 

Applying  the  classical  concept  of  energy  cascade,  (Tennekes  and  Lumley,  1983),  which 
assumes  energy  be  transported  from  larger  to  smaller  scales,  energy  transferred  from  grid 
scales  to  sub-grid  scales  could  considered  a  dissipation  of  energy;  the  task  of  a  SGS  model 
should  be  then  to  ensure  the  correct  ratio  of  dissipation  in  the  numerical  solution  of  the 
transport  equation  of  grid  ones  following  the  classical  energy  cascade  concept  (see  Fig.  2.3). 

Applying  the  classical  concept  of  energy  cascade,  Tennekes  and  Lumley  (1983) ,  which  as¬ 
sumes  energy  be  transported  from  larger  to  smaller  scales,  energy  transferred  from  grid  scales 
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to  sub-grid  scales  could  considered  a  dissipation  of  energy;  the  task  of  a  SGS  model  should 
be  then  to  ensure  the  correct  ratio  of  dissipation  in  the  numerical  solution  of  equation  2.8. 
In  this  frame,  the  classical  energy  cascade  concept  would  take  form  shown  in  Fig.  2.3. 


Mean  Flow 


Grid-Scale  Subgrid-scale 

Turbulence  Turbulence 

Fig.  2.3:  Kinetic  energy  exchange  between  grid-scales  and  sub-grid  scales 

The  simplest  way  to  model  a  dissipation  process  is  to  consider  it  analogous  to  molecular 
one  and  model  it  through  an  Eddy  viscosity  concept  (sub-grid  viscosity): 

nj  =  ^TkkSij  -  2 utSij  (2.18) 

It  is  to  be  noted  that  this  modeling  approach  assumes  the  turbulence  to  be  in  equilibrium, 
i.e.  production  of  turbulence  locally  balances  its  destruction.  In  other  words,  it  is  assumed 
the  small  scales  of  motion  adjust  instantaneously  to  the  perturbations  of  the  large  scales 
implying  the  existence  of  an  inertial  sub-range. 

2. 3. 1.1  Smagorinsky  Model 

By  far,  the  most  frequently  applied  SGS  model  is  the  one  developed  by  Smagorinsky 
(1963).  This  model,  based  on  dimensional  arguments,  assumes  an  Eddy  viscosity  based  on 
local  equilibrium.  The  model  is  based  on  the  following  quantities 


T{j  —  2  UtSij 

V,  =  C,A2|S| 
c  _ 

~  2  \aXj  dxj 

|S|  —  yJ^SijSij 


(2.19) 

(2.20) 

(2.21) 

(2.22) 


where  S{j  is  the  local  strain  rate  tensor  for  the  resolved  scales  velocity  field,  the  filter  size  A 
is  taken  to  be  proportional  to  the  local  grid  spacing,  and  Cs  is  an  empirical  proportionality 
constant,  which  depends  on  the  type  of  flow  being  modeled. 
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As  it  stands,  it  is  evident  that  the  Smagorinsky  model  suffers  from  two  drawbacks. 

First,  the  constant  C,  is  flow  and  grid  dependent  (Piomelli,  1997). 

Second,  the  model  does  not  account  for  the  rapid  decrease  of  length  scales  of  turbulence 
approaching  solid  surfaces,  and  thus  does  not  accurately  represent  the  physics  of  the  wall- 
layer.  Therefore,  a  correction  function  must  be  added  to  account  for  the  physical  decrease  of 
turbulence  level  nearing  a  solid  boundary. 

Wall  correction  usually  takes  the  form  of  a  van  Driest  damping  function  (eqn.  2.23). 


„«  =  CA2  (l  -  6-+/*+)3  |S|  (2.23) 

where 

z+  =  ^,  Ur  =lM’  =  25  (2.24) 

2. 3. 1.2  Structure  Function  Model 

Lesieur  (1987)  proposed  a  model  of  spectral  turbulent  viscosity.  The  idea  is  to  find  a 
relation  between  the  length  scale,  l,  and  the  cutoff  wave  number,  kc  =  (tt/A)  as  well  as  a 
relationship  between  the  velocity,  u,  with  the  cutoff  kinetic  energy,  E  {kc,  t).  A  transformation 
from  spectral  space  to  physical  space  is  achieved  using  the  resolved  second  order  Structure 
Function  for  the  velocity  vector  u,  which  can  be  expressed  as 

^7  =  (||ii(x,t)  -ii(a;  +  r,t)||2)  (2-25) 

which  could  be  related  to  local  resolved  strain  rate  and  vorticity  as 


§2  =  A2  (\S\  +  M2  ) 

1  /  dv,i  _  duj\ 

2  \  dxj  dxi  J 

\u\  =  yj&ij&ij 


The  corresponding  SGS  viscosity  can  be  expressed  as 

vSas  =  C„  A\/37  C„  =  O.lOSfc-3/2  (2-30) 


(2.26) 

(2.27) 

(2.28) 
(2.29) 
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Coefficient  CSf  is  obtained  assuming  a  Kolmogorov  cascade  and  is  considered  univer¬ 
sal  (Compte  and  Lesieur,  1998). 

The  resulting  Structure  function  (FS)  sub-grid  model  was  shown  (Compte  and  Lesieur, 
1998)  to  be  less  dissipative  than  the  Smagorinsky  model. 

An  important  improvement  of  this  concept  is  the  filtered  structure  function  (FSF)  (Compte 
and  Lesieur,  1998)  where  a  Laplacian  filter  T 


~  2 fij,k  +  )  + 

( fij+l,k  ~  2/ij,fc  +  fi,j-\,k)  + 
[fi,j,k+l  ~  2 fi,j,k  +  fij,k-l) 


(2.31) 


is  thrice  applied  to  the  velocity  field,  to  obtain  the  filtered  structure  function  Q ;  an^ 
corresponding  SGS  viscosity 


vsgs  —  Cfsf A 


(2.32) 


The  effect  of  the  filtering  is  to  enhance  the  contribution  of  the  smallest  resolved  scales, 
which  are  known  to  be  the  most  active  in  the  interchange  between  resolved  and  unresolved 
fields.  The  coefficient  Cfsf,  which  is  again  assumed  universal,  is  optimized  a  priori,  from 
benchmark  numerical  fields,  to  Cfsf  —  0.0084. 

Extensive  in-house  experience  on  FSF  model  (Hoffman,  1996)  has  shown  it  capable  to 
yield  quite  satisfactory  results  and  has  been,  therefore,  applied  in  the  present  study. 


2.3.1. 3  Dynamic  Model 

A  fundamental,  inherent  limitation  of  previous  discussed  SGS  models  lies  with  their 
dependence  on  a  coefficient  to  be  specified  a  priori,  on  the  basis  of  type  of  flow  and  grid 
(Smagorinsky)  or  a  priori  optimization  (structure  function).  This  drawback  could  be  over¬ 
come  applying  the  dynamic  procedure  proposed  by  Germano  et  al.  (1991)  which  allows  to 
compute  the  model  coefficient  (s)  dynamically  in  function  of  flow  field,  geometry  and  grid. 

Dynamic  procedure  was,  rightly,  heralded  as  a  major  advance  in  theory  and  understand¬ 
ing  SGS  modeling,  but  extensively  experience  on  different  versions  of  dynamic  Smagorinsky 
model  have  shown  that,  in  application  to  wall-bounded  flow  it  has  yielded  practical  results 
which  are  not  systematically  better  than  ones  yielded  by  basic  Smagorinsky  model  (Akselvoll 
and  Moin,  1995;  Balaras,  1995;  Cabot  et  al.,  1999;  Rodi  et  al.,  1997).  This  finding  was 
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confirmed  by  in-house  experience  (Laskowski,  1999;  Simons,  2000),  which  led  to  avoid  the 
application  of  this  procedure  for  the  present  investigation. 

2. 3. 1.4  Un-resolved  DNS 

It  was  often  remarked  in  the  literature  that,  in  case  of  2nd  order  Finite  Difference  dis¬ 
cretization  LES  equations,  the  contribution  of  SGS  has  the  dimension  of  the  leading  term  of 
the  discretization  error.  Furthermore  it  was  proposed  (Ghosal,  1995)  that  this  term  is  bigger 
than  the  one  of  model  SGS. 

If  this  is  the  case,  it  would  be  possible  to  perform  LES  without  SGS  model,  solving 
eqn.  2.8  with 


=  0  (2-33) 

and  leaving  the  numerical  error  provides  the  necessary  dissipation.  Under  the  numerical 
point  of  view  this  simulation  will  correspond  to  a  DNS  over  a  grid  too  coarse  to  resolve  the 
relevant  scales,  hence  the  naming  for  “un-resolved  DNS”  for  this  class  of  simulation. 

It  has  to  be  stated  that  practical  experience  of  present  authors  show  LES  to  be  sensible 
to  effect  of  the  SGS  model,  even  for  2nd  order  Finite  Differencing,  which  will  also  be  shown 
in  §  5.  On  the  other  hand,  it  is  actually  possible  to  apply  un-resolved  DNS  to  flows  close 
to  equilibrium  obtaining  realization  which  are  numerically  stable  and  comparable  with  full 
LES  (with  SGS  model).  Therefore  this  approach  represents  an  useful  benchmark  for  the 
sensibility  of  results  to  the  effect  of  SGS  model  and,  for  this  reason,  has  been  applied  in 
present  investigation. 

2.3.2  Self-similar  SGS  models 

In  is  well  known  (Hartel  and  Kleiser,  1998)  that  sub-grid  viscosity  cannot  fully  repro¬ 
duce  the  actual  behavior  of  sub-grid  scales  as  put  in  evidence  by  DNS  (Hartel  and  Kleiser, 
1994),  neither  the  effective  process  of  energy  transfer  between  the  resolved  and  unresolved 
scales.  This  process  can  better  understood  considering  the  transport  equation  for  the  resolved 
turbulent  kinetic  energy: 


dq2  dg2Uj  _  _d_ 
dt  dxj  dxj 


I  -2 puj  -  2rijUi  + 


J _df\ 

Re  dxj  J 


2  dui  dui  , 

Re  dxj  dxj  +  Tij  ij 


(2.34) 


where  the  last  two  terms  represent  the  production  and  destruction  of  resolved  kinetic  energy, 
whereas  the  other  terms  redistribute  energy  via  advection  and  diffusion  (Hartel  and  Kleiser, 
1994).  It  is  important  to  note  that  the  last  term  can  be  positive,  corresponding  to  a  reverse 
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flow  of  energy  from  the  small  scales  to  the  large  scales,  phenomenon  commonly  referred  to  as 
backscatter.  It  is  found  that,  in  case  of  wall-bounded  flows,  this  process  is  primarily  present 
and  important  in  the  buffer  region  of  the  wall  layer,  that  is  the  region  near  a  solid  surface 
where  7  <  z+  <  20.  Therefore,  sub-grid  viscosity  models  will  never  be  able  to  fully  reproduce 
the  behavior  of  sub-grid  scales  in  wall  layer  (Hartel  and  Kleiser,  1994). 

In  this  optic,  a  more  physical  model  for  interaction  of  grid  and  sub-grid  scales  could  be 
the  one  proposed  by  Hartel  and  Kleiser  (1994)  and  shown  in  Fig.  2.4. 


Grid- Scale 
Turbulence 


Subgrid-scale 

Turbulence 


Fig.  2.4:  Kinetic  energy  interchange  between  grid-scales  and  sub-grid  scales.  (P— turbulent 
production, E=turbulent  dissipation,  e= interchange) 


where  turbulent  kinetic  energy  can  flow  from  grid  scales  to  sub-grid  scales  (forward  scatter) 
and  from  sub-grid  scales  to  grid  scales  (backward  scatter). 

This  mechanism  could  be  taken  in  account  with  applying  self-similarity  concept  of  SGS 
model  (Meneveau  and  Katz,  2000). 

This  approach  is  based  on  the  observation  that  most  of  the  interchange  between  the 
resolved  and  unresolved  scales  happens  between  the  smallest  resolved  scales  and  the  biggest 
unresolved  scales.  Furthermore,  these  two  families  of  scales  exhibit  similar  behavior,  due  to 
the  fact  that  the  filter  which  artificially  separate  the  inertial  range  (Fig.  2.2)  and  is  the  actual 
size  of  the  filter  which  causes  a  scale  belong  to  one  of  these  two  fields. 

Self-similarity  of  smallest  of  the  resolved  scales  and  the  biggest  of  unresolved  scales  could 
be  in  evidence  with  second  filtering  operation  over  the  flow  field  (Fig.  2.5) 


(2.35) 


where  f  is  a  general  operator. 

Bardina’s  basic  proposal  (Bardina  et  al.,  1984)  is  to  filter  the  resolved  field 
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Fig.  2.5:  Region  of  interaction  between  smallest  resolved  scales  (1)  and  biggest  unresolved 
scales  (2) 


nj  =  Ui  Uj  -  Ui  uj  {2.Z6) 

This  approach  has  shown  its  capability  to  reproduce  the  physics  sub-grid  scales,  but  was 
found  unable  to  create  an  adequate  dissipation  at  the  level  of  grid  scales  (Mom  and  Mahesh, 
1998;  Piomelli,  1997),  and  usually  supplemented  by  dissipative  term  in  form  of  sul>grid 
viscosity,  to  obtain  the  “Mixed  model” 


Tij  =  Ui  Uj  —  U{  Uj  —  VSgs  Sij 


(2.37) 

(2.38) 


where  uSGs  can  be  obtained  with  Smagorinsky  model  or  through  dynamic  procedure  Moin 
and  Mahesh  (1998). 

Further  analysis  of  filtering  approach  led  the  conclusion  Domaradzki  and  Saiki  (1997), 
Domaradzki  and  C.  (1999);  Moin  and  Mahesh  (1998)  that  this  operation  should  be,  in  prin¬ 
ciple,  applied  to  the  fully  flow  field  m  because  main  interactions  are  localized  between  the 
smallest  resolved  scales  and  the  unresolved  scales  of  sizes  y  <  i  <  A. 

In  this  optics  Bardina’s  model  could  be  seen  to  approximate  total  Ui  on  the  basis  of  Ui 


Ui  «  Ui 


(2.39) 


Following  this  reasoning,  different  proposals  Domaradzki  and  Saiki  (1997);  Domaradzki 
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and  C.  (1999);  Shah  and  Ferziger  (1995)  have  been  made  to  obtain  an  estimation  ux  of  m 
from  ui . 

In  particular,  the  proposition  of  Shah  and  Ferziger  Shah  and  Ferziger  (1995)  is  to  apply 
to  a  top-hat  filter  to  ui  to  obtain  a  3-diagonal  systems  in  U{ 


(2.40) 


where  Ui  is  an  approximation  to  the  FULL  flow-field  U{. 

It  is  stated  (Shah  and  Ferziger,  1995)  that  this  formulation  produces  enough  dissipa¬ 
tion  to  make  unnecessary  the  SGS  diffusivity  term.  This  model  was  intensively  tested  in- 
house  Louda  (2001)  and  found  to  yield  better  representation  of  resolved  scales  than  the  one 
produced  by  sub-grid  viscosity  models  at  higher  Reynolds  number,  although  its  computa¬ 
tional  cost  makes  it  prohibitive  for  baseline  SGS  model  for  the  complete  study.  It  will  be 
applied  to  assess  the  effect  of  sub-grid  viscosity  model  on  simulation  of  coherent  structures. 
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3.  NUMERICAL  SOLUTION  OF  LES  EQUATIONS 
3.1  Time- Integration 

Paramount  influence  of  accuracy  of  discretization  on  LES  (Ghosal,  1995)  and,  in  con¬ 
sequence,  on  simulation  of  coherent  structures,  leads  to  the  need  of  discussing  in  depth  the 
numerical  formulations  applied  in  the  present  study. 

As  seen  in  2.1  ,  transport  equations  for  resolved  scales  take  the  form 


duj 

dt 


dui 

dxi 

=  0 

dr -(a)  _ 

(3.1) 

dUiUj 

dp  !  1  d2Ui 

(3.2) 

dxj 

dxi  Re  dxjdxj 

OXj 

to  be  solved  by  numerical  means. 

The  momentum  equation  can  be  re-written  in  functional  form 


dui 

dt 


-pL  +  Hfc.Ui) 


(3.3) 


where  H(ut, Uj)  contains  the  adjective  C(y.,,  u,)  and  the  diffusive  D(u j)  operators 


u, _ ,  3u , 

HiiH.Ui)  = —g— + 


1  &Ui 

dxj  V  R-e  9xj 


r-fj  =C(ui,  Uj)  +  D(ui) 


(3.4) 


In  the  present  code,  the  above  equations  are  solved  advancing  in  time,  applying  the 
second  order  accurate  Adams-Bashforth  algorithm  (Balaras,  1995;  Hoffman,  1996),  whic 
can  be  written  as 


u?+l  -Wj 

A  t 


=  ^H(uin,ui 


n  -  ij m’-wr1) 


dxi 


(3.5) 


This  formulation  involves  four  unknowns  at  the  time  level  n+1:  u,v,w,  and  p,  and 
only  three  transport  equations  for  the  three  components  for  momentum.  Matter  of  fact, 
in  incompressible  flows  the  unsteady  Navier-Stokes  equations  involves  temporal  derivatives 
in  the  momentum  equation  but  not  in  the  mass  equation.  In  consequence,  the  divergence 
of  the  velocity  field  appears  only  as  a  constraint  on  the  solution  of  the  momentum  equa¬ 
tions.  This  constraint  is  satisfied  during  time  integration  applying  the  fractional  time  step 
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approach  (Chorin,  1968;  Ferziger  and  Peric,  1997).  This  procedure  consist  of  two  successive 
steps:  predictor  and  corrector  step. 

To  estimate  the  velocity  field  at  time-step  n  +  1,  the  predictor  step  computes  an  inter¬ 
mediate  velocity  field  Uin+*,  applying  the  pressure  field  at  time  level  n 


tt  n+frac\2  _  rr  n  o  1  ,  , 

± — sr3-  =  V)  -  -2H{uT\ url) 


gr 

dxi 


(3.6) 


the  resulting  uT+//rac12  is  not  divergence-free. 

The  final,  divergence  free,  solution  is  imposed  in  the  corrector  step,  which  is  the  difference 
of  the  two  previous  equations: 


u. 


A  t 


_g_ 

dxi 


p")  =  - 


dsp 

dxi 


(3.7) 


while  the  divergence-free  constraint  is  applied  through  a  pressure  Poisson  equation  whose  the 
source  term  is  supplied  by  the  divergence  operator  for  the  intermediate  velocity 


V  •  izr+1  =  V  •  izr+*  ~  AfV  •  =  0  (3-8) 

Imposing  the  final  solution  to  be  divergence-free: 

V  •  uf +1  =  0  (3.9) 

the  above  can  rewritten  as 

;  t 

d25p  1  du?+$  (3 10) 

dxi 2  A t  dxi 

to  find  U  and  solve  the  corrector  step. 

Finally,  the  pressure  field  is  updated 

p"+l  =pn  +  Sp  (3.11) 

In  this  form  the  velocity  is  explicitly  solved  whereas  the  pressure  is  solved  implicitly. 
It  should  again  be  noted  that,  in  this  frame, p  is  no  more  the  adimensional  form  of  the 
thermodynamic  pressure  ,  but  a  modified  pressure  containing  the  isotropic  part  of  the  SGS 
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stress  tensor.  Solving  this  modified  pressure  implicitly  does  not  require  the  knowledge  of  the 
isotropic  part  of  the  SGS  tensor  but  in  the  same  time  make  more  difficult  compare  results 
based  on  the  pressure  field.  As  the  velocity  is  solved  through  a  fully  explicit  algorithm  (for 
viscous  and  advective  terms),  the  permissible  time  step  is  limited  by  stability  considerations 
based  on  both  convection  (C)  and  diffusion  (D).  Typical  limits  for  At  for  Adams-Bashforth 
algorithm  for  present  application  are 


Convection 

Diffusion 


A  te 


<  0.35 


<0.1 


(3.12) 

(3.13) 


where  VT  is  the  sum  of  the  molecular  and  SGS  viscosity  (obviously  the  SGS  viscosity  goes  to 
zero  for  DNS  simulations). 


3.2  Spatial  Discretization 

In  present  investigation,  LES  equations  are  solved  on  stretched,  orthogonal  grids.  Stag¬ 
gered  grid  is  adopted  to  avoid  the  generation  of  spurious  pressure  oscillations  associated  with 
the  solution  of  the  incompressible  Navier-Stokes  equations  on  collocated  grids  (Ferziger  and 
Peric,  1997).  Fig.  3.1  shows  the  locations  of  u,  v,  wandp:  pressure  is  located  at  the  cell  center 
and  the  velocity  components  are  located  at  the  corresponding  flux  surfaces. 


vr(i,j,k) 

v(i,j,k) 

* 


_U  (i,jrk) 


Fig.  3.1:  Location  of  the  velocity  coordinates  and  pressure  on  a  staggered  grid 

The  Arakawa  (1966)  form  of  the  advection  term  is  adopted  to  minimize  aliasing  errors, 
which  result  from  misrepresentation  of  shorter  waves  by  the  grid  (Hoffman,  1996). 
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(3.14) 


dul  Uj 
dxj 


1  fduiuj  duj  g*7 
dxj  J  dxj  'dxj 


This  formulation  imposes  an  integral  constraint  of  mean  kinetic  energy  by  treating  the 
advective  terms  as  the  average  between  the  divergence  and  conservative  form.  In  the  most 
common  approach,  spatial  derivatives  are  discretized  applying  2nd  order  central  finite  differ¬ 
ences.  Staggered  set-up  needs  the  definition  of  velocities  on  the  faces  of  cell  velocity  as  shown 
in  Fig.  3.2. 


Ax . 
i 


u-velocity 
control  volume 

v-velocity 
control  volume 


pressure 
control  volume 


Fig.  3.2:  Cell  and  flux  surfaces  for  velocity  component  u{ 

The  advective  term  is  discretized  as 

dUjUj  _  W|j~  Uj\t  ~  ui\j  Uj\j  /g  JCj) 

dxj  A  Xj 

where  u]\f  denotes  the  interpolation  of  the  j  component  of  velocity  in  the  i  direction  at  the 
±  flux  surface.  For  the  present  investigation,  linear  interpolations  are  applied,  to  maintain 
2nd  order  accuracy  for  the  derivative  at  the  velocity  location. 

Following  the  same  approach,  the  molecular  diffusive  fluxes  are  discretized  as 

d2uj  _  d  I  m 

dxjdxj  dxi  ^  dxj 

where  2nd  second  order  discretization  is  used  for  the  first  derivative  ^  at  flux  surfaces. 


dul 

dxj 


(3.16) 
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Finally,  the  pressure  gradient  is  discretized  as: 


dp  _  Pi+l,j,k  Pi,j,k 

dXi  i+Lj'k  Xi+U,k  ~  Xhi’k 


(3.17) 


The  limitations  of  second  order  discretization  for  LES  are  well  known  (Baggett  et  al., 
1997).  Therefore  a  fourth  order  discretization  was  also  applied  to  assess  the  influence  of 
numerical  accuracy  on  simulation  of  coherent  structures. 

In  present  approach,  fourth  order  accuracy  is  obtained  applying  the  procedure  proposed 
by  Verstappen  and  Veldman  (1997),  where  advective  and  diffusive  flux  are  computed  no  only 
on  the  basic  cell  (Fig.  3.3),  but  also  a  superposed  one,  thrice  bigger,  one. 


- 

PI 

fit  M 

Fig.  3.3:  Small  and  large  grid  for  4t/!  order  discretization 

Linear  interpolations  are  applied  on  both  grids,  corresponding  to  2"d  order  accuracy  on 
both  grids. 

SJlf  =  0.5(-  +  t75)  (3-18) 

Summing  up  the  Right  Sides  calculated  on  the  two  grids,  the  final  form  is  obtained 

?*  =  *(-&.  + HI*,*)}  +  (319) 

dt  8  V  dxi 


smallgrid 


4(-S+^’5i) 


widegrid 


which  yields,  on  weak  stretched  grids,  fourth  order  accuracy  for  final  solution. 

The  method  proposed  (Verstappen  and  Veldman,  1997)  is  based  on  the  Richardson 
Extrapolation.  Considering  a  second  order  central  difference  derivative  on  the  usual  stencil 

of  3  nodes 


A*,-)  =  /(— + 


(3.20) 
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and  a  generic  stencil  of  2  *  j  +  1  nodes 


J &->).  + (3.21) 

the  Richarson  extrapolation  leads 

j_  1^2+AzMBzA _ 1  ffa+i)  ~  +  q2  (a*)  (3.22) 

f{  l)  f-1  2AX  j2- 1  2jAx  {  x) 

Spectral  behavior  of  the  new  formulation  depends  on  the  size  of  the  new  grid,  its  per¬ 
formances  rapidly  decaying  with  the  increasing  dimensions  of  the  new  grid.  Remembering 
the  implicit  filtering  performed  from  the  finite  differencing,  eqn.  3.22  is  equivalent  to  a  new 
transfer  function  G(r )  (Fig.  3.4) 

,  9f{xi+l)  1  —  /(xi-s) 

f[Xi)  "8  2AX  8  6AX 

f{xi)  =  G{y-Xi)f{y)dx 

=  4^(A--|r|)  +  l3(2bJ?(3A'-|r|) 

used  to  compute  the  filtered  energy  spectra 

Wi{k)  =  \G"(k)\2En(k)  (3.26) 

where  G"  (k)  is  the  Fourier  pair  of  G(r). 

This  form  of  discretization  has  been  intensively  tested  in  house  tests  and  improved 
simulation  of  resolved  turbulence  was  found  for  the  lower  range  of  Re  examined.  At  higher 
Re  the  velocity  spectra  present  a  deviation  from  the  normal  behavior  that  prevents  its  use. 

Both  2nd  and  4t/l  central  schemes  do  not  satisfy  the  monotonic  property  (see  Hirsh  (1990) 
for  a  classical  discussion  of  the  matter)  and  may  the  introduce  spurious  minima  and  maxima 
(“wiggles”)  in  the  solution  when  applied  to  flows  at  at  high  Re.  The  phenomenon  can  be 
minimized,  or  avoided,  implementing  upwind  or  upwind-biased  discretizations  (Hirsh,  1990) 
where  the  computational  pencil  is  weighted  in  upstream  direction.  In  the  application,  an 
upwind-biased  3rd  order  (Kawamura  et  al.,  1984)  scheme  is  introduced.  This  scheme  is  made 
of  the  sum  of  a  5  points  4th  order  central  discretization  CAth(fi- 2,  fi-i,  fi,  fi+ 1)  fi+ 2)  and  a  4 
order  dissipation  D4th(fi- 2,  /*-i,  /*,  fi+i,  fi+ 2) 

~  =  Cm(fi- 2,  fi- 1,  fi,  fi+ufi+2)  +  ctD4th(fi-2,  fi-u  fi+ufi+2)  (3-27) 


(3.23) 

(3.24) 

(3.25) 
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Transfer  Function  G" 


Fig.  3.4:  Spectral  behavior  of  G(r) 


where  0  <  a  <  1  with  a  =  0  yields  central  4th  order  scheme  and  a  =  1  yields  the  full  3rd 
order  up-wind  scheme.  While  this  scheme  is  not  strictly  monotonic,  it  was  found  adequate  to 
control  the  formation  of  wiggles  in  present  applications  Hoffman  (1996);  Giammanco  et  al. 
(2001). 

3.3  Multi-domain  technique 

Application  of  LES  to  flow  around  bodies,  such  as  cylinder  to  be  considered  for  the 
present  study,  poses  challenging  fundamental  problems  concerning  the  implementation  of 
wall  condition  and  economical  solution  of  equation  for  the  pressure. 

In  fact,  Fast  Poisson  Solvers  (FPS)  (Schmidt  et  al.,  1984)  cannot  directly  be  applied 
to  this  case  and  capacitance  matrix  technique,  currently  applied  to  extend  FPS  to  complex 
geometries  (Schmidt  et  al.,  1984;  Akselvoll  and  Moin,  1995),  requires  the  solution  of  two 
Poisson  problems  over  an  enlarged  domain  including  the  actual  computational  field. 

To  overcome  this  problem,  in  present  study,  LES  equations  are  solved  applying  the  multi- 
domain  (MD)  technique.  In  MD  approach,  the  computational  field  is  sub-divided  in  a  set 
of  non-overlapping  sub-domains  and  the  original  problem  is  decoupled  in  a  set  of  two  sub¬ 
problems:  a  sub-domain  problem  (consisting  of  a  sub-set  of  decoupled  problems,  Fig.  3.5) 
and  an  interface  condition  (which  ensures  the  appropriate  coupling,  Fig.  3.6). 

The  advantages  of  this  procedure  lie  with  the  easiness  of  the  description  and  discretiza¬ 
tion  of  complex  geometries,  good  grid  control,  feasibility  of  parallel  execution  on  multi-node 
computers  (Simons,  2000),  together  with  the  possibility  to  work  on  simple  sub-domains  where 
Fast  Poisson  Solvers  (FPS)  could  be  applied. 
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1 


Fig.  3.5:  Example  of  subdivision  in  3  sub-domains  f2i,fi2,n3  connected  by  interfaces  ri,r2 


Fig.  3.6:  Boundary  points  update  at  vertical  interface  between  sub-domains  Qi,Qj 


The  discretization  of  pressure  in  a  multi-domain  formulation  is  straightforward;  it  is 
nevertheless  necessary  to  introduce  some  nomenclature:  letxn  and  xr  denote  the  solution 
vectors  corresponding  to  the  sub-domains  and  to  the  interfaces  whose  collection  adds  up  to 
the  original  domain.  In  a  collocated  formulation  the  interface  T  is  not  well  defined  at  discrete 
level  (viz.  there  are  no  unknowns  on  the  boundaries  separating  two  sub-domains).  Here  T 
is  defined  as  one  (or  more)  layers  of  unknowns  next  to  the  boundary  which  are  therefore 
eliminated  from  sub-domain  count,  where:  x  =  xnUxr- 

Fourier  expansion  has  been  used  for  unknowns  and  right  hand  sides  of  the  equations 
involved  in  the  homogeneous  span-wise  direction  y,  so  that  the  original  three-dimensional 
Poisson  problem  is  reduced  to  a  set  Ni  of  two  dimensional  Helmholtz  problems,  which  differ 

from  each  other  exclusively  for  a  parameter  Xm,  which  essentially  depends  on  m  =  1,  Ni. 

;  2 

It  must  be  noted  that,  because  of  the  symmetry  of  the  Fourier  transformed  equations, 
Ni  problems  need  to  be  solved.  Ordering  the  unknowns  in  a  domain  wise  fashion,  followed 
by  the  inter-facial  unknowns,  and  applying  block  Gaussian  elimination  to  the  original  system 
of  equations  arising  from  the  discretization,  brings  to  the  following  system 
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(  014  ...  0  ai,iv  •••  ai,jv+iw  ^ 

(  xni  } 

bi  \ 

0  ...  Ow.JV  aJV,W+l  •••  <3iV,jV+M 

X&N 

bN 

0  0  0  sn+i,n+i  •••  SN+1,N+M 

Xri 

bN+  i 

0  0  0  . 

•• 

•• 

^  0  0  0  sn+m,n+i  -  sn+m,n+m  / 

V  XrM  / 

\  bi v+M  / 

(3.28) 


where  aiti  is  the  matrix  representing  the  discretization  of  the  left  sides  of  in  N  and  M 
are,  respectively,  the  number  of  sub-domains  and  interfaces  number  of  unknowns,  and 
bij  are  defined  as 


N 


Sij  =  Aij  -  aikak£akj  Vi,  j  =  N  +  1,  N  +  M 

k- 1 

(3.29) 

N 

bj  =  bij  -  ^2 ajkakkbk  Vi,  j  =  N  +  1,  N  +  M 

(3.30) 

k—l 

(3.31) 

The  interface  unknowns  can  be  determined  from  the  lower  right  sub-block  of 


S  Xy  = 


S/V+l.JV+l 
Sn+M,N+  1 


Sn+\,N+M 

Sn+M,N+M 


XTl 


X?M 


(3.32) 


The  sub-matrix  s  is  generally  referred  to  as  the  Schur  complement  matrix. 

The  connection  between  the  sub-domain  and  interface  unknowns,  which  is  apparently 
lost,  is  now  embedded  in  the  right  side  of  equation  above. 

To  summarize,  the  solution  algorithm  reduces  to  the  following  a  three-step  procedure: 


1.  Separate  solution  of  an  elliptic  problem  in  each  of  the  sub-domains  Q*,  i  =  1,  iV,  using 
the  original  boundary  conditions  at  the  external  boundaries,  homogeneous  Dirichlet 
boundary  conditions  at  the  internal  boundaries  and  its  original  right  sides. 

2.  Assembling  of  the  right  sides  of  the  inter-facial  problem  using  the  previous  step  values, 
and  solution  of  the  Schur  complement  problem. 

3.  Separate  solution  of  an  additional  elliptic  problem  in  each  of  the  sub-domains  f2*, using 
the  original  boundary  conditions  at  the  external  boundaries,  original  right  sides  6*. 
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Although  the  overall  cost  of  the  elliptic  kernel  appears,  at  a  first  sight,  more  than  doubled 
(two  elliptic  problems  plus  an  interface  problem  of  smaller  size),  this  overhead  is  largely 
compensated  by  the  fact  that  this  approach  not  only  allows  to  recover  the  use  of  fast  Poisson 
solvers  on  rectangular  domains  (to  tackle  first  and  last  step),  but  also  opens  the  way  to  parallel 
execution. 

Additionally,  in  an  iterative  solution  context,  it  provides  the  opportunity  of  relaxing  the 
stiffness  of  the  full  scale  original  problem  (which  is  reduced  to  a  set  of  smaller  size,  better 
conditioned  sub-problems)  and  is  intrinsically  suited  for  zone  approaches  and  heterogeneous 
problems. 

For  application  of  the  LES,  which  effectively  requires  the  solution  of  several  elliptic 
equations  with  varying  right  sides  for  a  large  number  of  times  steps,  it  can  be  argued  that 
it  is  convenient  and  feasible  to  compute  and  store  the  inverse  of  the  Schur  matrix  through  a 
pivoted  LU  factorization  procedure  in  a  pre-processing  step. 

The  overhead  of  the  factorization  will  be  quickly  amortized  over  the  large  number  of 
time  steps  necessary  to  build  the  statistics.  It  has  to  be  underlined  that,  when  the  original 
problem  is  rank  deficient,  it  will  show  up  in  step  2  of  the  MD  algorithm  and  has  therefore  to 
be  tackled.  Present  code  keeps  the  capability  to  solve  the  discretized  elliptic  problems  found 
in  first  and  last  step  in  a  direct  or  iterative  fashion  (Simons,  2000). 

Extensive  previous  testing  on  simulation  of  boundary  layer  (Schmidt  et  al.,  1984)  or 
flow  over  a  backward  looking  step  (Simons  et  al.,  2001)  has  put  in  evidence  the  efficiency  and 
accuracy  of  the  procedure  and  it  yielded  satisfactory  results  in  the  present  study  (Giammanco 
et  al.,  2001). 
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4.  COHERENT  STRUCTURE  CONCEPT 


4.1  introduction 

Introduction  of  coherent  structure  concept  in  the  study  of  physics  and  phenomenology 
of  turbulent  flows  was  led  by  the  observation  (Brown  and  Roshko,  1975)  that  the  ensemble 
of  instantaneous,  three-dimensional  motions,  which  make  turbulence  (Tennekes  and  Lumley, 
1983),  contains  events  larger  and  longer  lived  than  the  scales  of  turbulent  transport  (Lesieur, 
1987;  Tennekes  and  Lumley,  1983). 

These  events  were  shown  to  have  a  level  of  internal  organization  much  higher  than 
turbulent  background,  leading  to  concept  of  coherent  structure.  Even  if  different  definitions 
of  coherent  structure  can  be  found  in  literature,  the  one  proposed  by  Hussain  (1986)  and 
already  cited  in  chapterl  1  was  found  adequate  for  present  purposes 

A  coherent  structure  is  a  connected  turbulent  fluid  mass  with  instantaneously 
phase  correlated  vorticity  over  its  spatial  extend.  . . .  We  choose  to  designate  the 
instantaneously  space  and  phase-correlated  vorticity  as  coherent  vorticity.  Thus, 
coherent  vorticity  is  the  primary  identifier  of  coherent  structures,  which  have  dis¬ 
tinct  boundaries  and  independent  territories. 


As  implied  by  the  above  definition,  vorticity  is  considered  the  main  physical  parameter 
characterizing  the  phenomenon  (Farge  et  al.,  1996). 

Paramount  importance  of  coherent  structures  lies  with  the  growing  understanding  of 
the  existence  of  an  underlying  order  in  chaos-looking  turbulence,  the  difficulty  lying  with 
identifying  and  assessing  of  respective  effects  of  deterministic  organization  and  chaos-like 
background  and  interaction. 

...frustration  results  largely  from  the  mixture  of  chaos  and  order  and  the  wide  range  of 
length  and  time  scales  that  turbulent  flows  possess.  (Moin  and  Mahesh,  1998). 

Coherent  structures  or  quasi- coherent  (this  being  the  case  for  organization  observed  in 
space  only  (Robinson  et  al.,  1990))  structures  have  been  identified  in  a  variety  of  turbulent 
shear  layers.  This  type  of  organization  is  more  evident  in  free  layers,  where  most  of  flow 
activity  re-organizes  itself  in  very  large,  quasi-periodic,  vortices  such  as  those  can  be  observed 
for  mixing  layers  (Brown  and  Roshko,  1975),  plane  wakes  (Hussain  and  Hayakawa,  1987), 
wakes  behind  a  solid  body  (Lyn  et  al.,  1995)  and  can  be  put  in  evidence  by  numerical 
simulation  (Fig.  4.1). 

Concerning  wall  layers,  a  variety  of  quasi-coherent  structures  was  put  in  evidence  for  the 
case  of  attached  boundary  layer  (Robinson  et  al.,  1990;  Spalart,  1988),  and,  in  particular, 
the  well  known  horse-shoe  or  hair-pin  (depending  on  Reynolds  number  of  flow)  vortex  has 
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Fig.  4.1:  Vortex  structures  in  wake  behind  cylinder  of  square  cross-section,  Re/r=22,000. 
LES 

drawn  attention  of  researchers,  being  considered  the  typical  structure  associated  to  turbulence 
production  in  buffer  layer  (Fig.  4.2) 

However,  recent  DNS(s)  tend  to  indicate  that  the  preferred  organization  in  buffer  and 
inertial  layer  consists  of  structures  elongated  in  stream-wise  direction  (Chong  et  al.,  1998; 
Jeong  et  al.,  1997).  This  conclusion  seem  to  be  confirmed  by  recent  in-house  DNS  (Rambaud, 
2000)  for  plane  channel  at  equilibrium,  where  “true”  horse-shoe  vortex  seem  to  be  relatively 
scarce  (Fig.  4.3) 

In  case  of  separate  wall  flows,  quasi-coherent  structures  have  been  put  in  evidence  in 
case  of  backward  facing  step  (Scarano  et  al.,  1999). 

Overall,  the  need  to  understand  the  physics  of  coherent  structures  and  their  influence 
on  the  average  behavior  of  the  flow  is  evident  for  all  classes  of  flows. 

To  quote  two  opposite  views  of  the  matter:  Hussain  (1986)  concludes  that,  in  case  of 
developed  turbulence,  behavior  of  coherent  and  un-coherent  fields  is  similar,  while  Farge  et  al. 
(1996)  takes  the  opposite  view  and  proposes  that  understanding  and,  possibly,  modeling  the 
different  physical  behavior  of  coherent  structures  and  un-coherent  background  might  be  the 
key  to  understanding  and  successful  prediction  of  turbulence.  In  particular,  Farge  considers, 
on  the  basis  of  behavior  of  two-dimensional  turbulence,  that  turbulent  flow  could  be  sub- 


Fig.  4.2:  Horse-shoe  structure  detected  in  plane  channel  at  ReT— 180.  DNS  (Rambaud,  2001) 

divided  in  vorticity- dominated  regions,  strongly  sheared  regions  and  weakly  sheared  regions, 
each  of  them  characterized  by  different  energy  spectra,  whose  different  characteristics  are 
masked  by  statistical  approach  embodied  by  classical  Kolmogorov  spectrum  (Tennekes  and 
Lumley,  1983). 

The  considerations  lead  to  the  need  to  understand,  not  only  shape  and  dynamics  of 
typical  structure,  but  also  the  quantitative  effect  of  coherent  structures  on  the  overall  behavior 
of  turbulent  field. 

This  being  the  case,  it  would  appear  that  the  numerical  simulation  would  be,  the  most 
appropriate  approach  to  this  task,  because  it  offers  the  possibility  to  access,  identify  and 
study  all  structures  which  might  be  present  within  flow-field  under  study. 

4.2  Deduction  of  coherent  structures 

4.2.1  Detection  criterion 

Following  the  arguments  presented  in  the  previous  paragraph,  the  basic  problem  would 
be  to  detect  organized  structures  within  three-dimensional  instantaneous  flow-fields  produced 
by  DNS  and  LES.  It  is  clear,  from  the  previous  discussion,  that  any  possible  detection  criterion 
would  involve  the  gradients  of  velocity  vector  Therefore,  the  starting  point  should  be  the 
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Fig.  4.3:  Vortex  structures  detected  in  plane  channel  at  Rer  =  180;  Top  view.  DNS  (Ram- 
baud,  2000). 

gradient  tensor  A  of  u*. 

Applying  analysis  and  nomenclature  proposed  by  Cantwell  (1992)  for  A  the  basic  oper¬ 
ator  to  be  considered  would  be 

A  =  Viii ;  Aij  =  (4.1) 

which  can,  obviously,  be  decomposed  in  symmetric  rate  strain  tensor  S 

S{j  =  —  (A{j  +  Aji)  (4.2) 

and  anti-symmetric  rate  of  rotation  Q 

^ ij  ~  2  —  Aj*)  (4.3) 

The  eigenvalues  A  of  A  are  solution  of  the  characteristic  equation 
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A3  +  PA2 +  QA  +  X  =  0 


(4.4) 


where 

P  =  -Su  (4-5) 

Q  =  (4.6) 

R  =  I(-P3  +  3PQ  -  SijSjkSki  -  (4.7) 

o 

Within  this  frame,  three  different  detection  algorithms  can  be  found  in  present  literature 

1.  “D”  (Cantwell,  1992)  criterion  :  in  this  case  the  discriminant  D 

D  =  ^ R 2  +  Q 3  (4.8) 

t: 

of  the  above  characteristic  equation  is  tested.  The  region  where  D  is  strictly  positive 
is  supposed  to  contain  the  coherent  structures. 

2.  “Q”  (Hunt  et  al.,  1988)  criterion:  in  this  case,  the  region  where  the  second  invariant  of 
Viii 


g  =  i/2(||5||2-||fi||2)1 


(4.9) 


is  positive. 

3.  A2  criterion  (Jeong  et  al.,  1997):  in  this  case,  a  connected  region  where  two  of  the  three 
eigenvalues  A  of  the  operator  S 2  +  fi2  are  negative  is  defined  as  “vortex  core”;  assuming 
Ai  >  A2  >  A3,  this  definition  leads  to  the  criterion 

A2  <  0  (4.10) 

for  the  detection  of  a  “vortex  core”.  Focusing  on  the  “vortex  core”  allows  to  separate 
true  coherent  vortices  from  other  possible  (incoherent)  regions  of  high  Q. 

It  was  pointed  out  by  a  very  recent  analysis  (Rambaud,  2000),  the  three  previous  criteria 
correspond  to  three  different  concepts  of  organized  vorticity 

•  D  detects  the  regions  where  stream-lines  are  helicoidal 
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•  Q  detects  the  region  where  local  vorticity  are  stronger  than  local  shear 

•  A2  detects  the  regions  of  minimum  pressure  within  the  body  of  structure. 


As  found  in  (Rambaud,  2000)  and  in  present  study,  these  criteria  do  not  necessarily  yield 
coinciding  results,  as  it  will  be  discussed  in  Chapter  5. 

4.3  Detection  procedure 

Under  theoretical  point,  all  the  above-introduced  criteria  take  the  form  of  a  test  on  the 
sign  (positive  or  negative,  depending  on  particular  criterion  selected)  of  a  testing  quantity, 
function  of  velocity  gradients  of  the  flow.  This  theoretical  criterion  must  be  integrated  into 
a  general  procedure  for  detection  and  definition  of  organized  structures  in  three-dimensional 
numerical  flow-fields.  Development  and  inner  working  of  the  resulting  algorithm  are  presented 
and  exemplified  for  the  case  of  LES  of  plane  channel  at  equilibrium,  at  Rer  =  180,  discretized 
over  a  grid  of  96x128x64  points  (Annex  A);  Q  criterion  is  applied  for  purpose  of  presentation, 
but  present  processing  and  arguments  lead  to  similar  results  applying  A2  criterion,  while  D 
criterion  yields  remarkable  different  results  ,  as  it  will  be  shown  in  Chapter  5. 

The  actual  construction  of  structures  in  the  simulated  flow-field  depends  of  discretization 
employed  for  the  resolution  the  flow-field.  In  the  case  of  present  code  (chapter  3),  this  is 
performed  by  finite  differencing  over  staggered  grid  where  the  three  components  of  velocity 
and  the  pressure  are  evaluated  at  different  positions  in  the  space  (Fig.  3.1  in  3.2). 

Since  the  pressure  node  is  located  at  the  center  of  the  computational  cell,  it  is  natural  to 
apply  the  detection  algorithm  at  the  pressure  node  and  evaluate  at  this  location  the  velocity 
derivatives  used  by  the  detection  algorithm.  Therefore,  the  initial  result  is  a  three-dimensional 
matrix  Q(x,y,z )  located  at  the  pressure  nodes. 

Theoretically,  all  the  grid  points  where  Q  >  0  should  be  considered  part  of  a  coherent 
structure,  but,  under  the  practical  point  view,  strict  application  of  this  criterion  yields  noisy 
fields,  where  it  is  impossible  resolve  the  individual  structures,  as  it  can  be  seen  in  Fig.  4.4. 

In  order  to  avoid  numerical  noise  created  by  values  close  to  0  it  is  necessary,  as  attested 
by  current  literature  (Baggett  et  al.,  1997),  to  introduce  a  threshold  value  (“trigger  level”, 
Tr),  sufficiently  higher  than  0,  above  which  the  criterion  is  considered  satisfied. 

Choice  of  a  trigger  level,  obviously,  retains  a  degree  of  arbitrariness,  but  a  sound  criterion 
could  be  deducted  from  the  analysis  of  the  probability  function  (PDF)  of  Q  (Fig.  4.5). 

It  can  be  observed  from  Fig.  4.5  that  positive  values  of  Q  for  global  PDF  span  from  0  to 
4000  (in  wall  units),  while  the  region  Q  <  10  can  be  interpreted  as  a  region  of  “white  noise” 
which  is  strongly  affected  by  precision  of  the  simulation.  Therefore,  an  useful  trigger  level 
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Fig.  4.4:  Iso-surfaces  of  Q  =  5;  plane  channel  at  equilibrium,  at  Rer  =  180;  lower  half 
channel. 
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Fig.  4.5:  PDF  of  Q  for  plane  channel  at  equilibrium,  at  ReT  =  180.  Zone  1:  0  <  z+  <  30; 
Zone  1:  30  <  z+  <  100;  Zone  3:  *+  >  100 


should  be  chosen  above  this  range.  As  title  of  example,  the  effect  of  trigger  levels  Tr  =  91 
(Fig.  4.6)  and  (Fig.  4.7)  Tr  =  200  are  presented  and  it  can  be  remarked  that  this  choice  has 
paramount  effect  on  the  final  result 

It  can  be  remarked  that  introduction  of  a  trigger  produces  distinct  and  separate  struc¬ 
tures  which  could  be  identified  and  defined  by  an  appropriate  algorithm. 

Further,  the  drastic  decrease  of  population  of  structures  with  increasing  Tr  is  evident. 
This  effect  can  be  appreciated  applying  the  statistical  tools  to  be  presented  in  next  paragraph: 
the  percentage  of  volume  of  computational  field  where  Q  >  Tr  for  the  two  cases  is  summarized 
in  Tab.  4.1,  where  a  dramatic  change  of  detected  structures  is  evident. 
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Fig.  4.6:  Iso-surfaces  of  Tr  =  91,  plane  channel  at  equilibrium,  at  Rer 
channel. 


180;  lower  half 


Fig.  4.7:  Iso-surfaces  of  Tr  =  200,  plane  channel  at  equilibrium,  at  ReT  =  180;  lower  half 
channel. 

Table  4.1:  Plane  channel  at  equilibrium  at  ReT  =  180  (Example).  Nominal  structures 

percentage  in  volume  _ 

Trigger  level  Tr  Volume  Q  >Tr 
91  1L33 

200  5.1376 


It  has  to  be  remarked  that  values  quoted  in  Tab.  4.1  correspond  to  the  volume  having 
Q  >Tr  without  discerning  whether  this  volume  is  made  of  only  true  structures  or  numerical 
noise  as  well;  this  matter  will  be  discussed  in  next  section. 

The  conclusion  to  be  drawn  is:  selecting  the  lower  value  of  trigger  level,  it  is  possible  to 
find  most  of  the  structures  possibly  present  in  the  field,  while,  applying  the  higher  level  only 


the  ones  corresponding  to  strongest  activities  turbulent  motions  are  retained. 

In  the  second  case,  a  sub-set  of  vortex  structures  is  identified,  whose  characteristics  are 
different  of  the  overall  average.  Structures  of  this  last  family  are  mostly  located  in  inner 
wall-layer  (z+  <  50),  as  it  can  be  seen  comparing  the  distribution  of  structures  across  the 
height  of  the  channel  for  Tr  =  91  and  Tr  =  200  (see  Fig.  4.8). 

It  is  important  at  this  point,  to  stress  the  following  concept:  all  the  following  plots, 
which  show  statistical  properties  such  as  single  PDF,  joint  PDF  are  based  on  histograms.  An 
histogram  can  be  defined  in  different  ways  but  the  simplest  one  consists  in  choosing  a  series 
of  equal-spaced  data  defining  an  interval  called  sampling  bin.  Considering  a  fixed  interval  to 
cover,  higher  is  the  number  of  points  chosen  to  represent  it,  smaller  is  the  size  of  the  bin  and 
smaller  is  the  number  of  points  that  fell  into  this  interval. 

Using  an  histogram  is  nothing  else  than  smoothing  the  raw  data;  this  concept  must  be 
kept  present,  when  examining,  for  example,  at  Fig.  4.8,  where  the  number  of  points  chosen 
and  the  interval  to  cover  are  clearly  mentioned.  The  bin  size  is  easily  extrapolated  from  these 
two  data. 

At  title  of  example,  in  aforementioned  figure,  ordinate  axis  represents  a  percentage.  It 
represents,  in  detail,  the  percentage  of  structures  which,  in  average,  have  their  barycenter 
within  the  limits  identified  by  the  corresponding  abscissa  and  the  bin  width.  Decreasing  the 
number  of  points  would  have  meant  increasing  the  bin  width  and  therefore  increasing  the 
percentage  associated  to  the  same  abscissa. 

In  any  case,  the  discrete  sum  of  the  ordinates  of  the  nodes  it  will  always  be  100,  unless 
the  plot  is  showing  a  zoom  of  a  region  of  interest,  which  is  usually  the  case  for  plot  of  vorticity, 
enstrophy  and  other  dynamic  quantities. 

In  the  present  investigation,  the  number  of  points  has  been  chosen  high  enough  to  mimic 
the  real  behavior  of  the  data,  but  with  a  bin  size  high  enough  to  perform  some  smoothing  of 
the  results. 

In  the  following,  the  number  of  bins  taken  and  the  interval  they  represent  will  be  quoted 
in  all  plots,  except  when  such  information  is  not  necessary. 

Coming  back  to  the  examination  of  the  two  families  of  structures  defined  by  the  two 
different  trigger  levels,  it  can  be  observed  that  they  have  different  physical  properties,  which 
is  shown,  for  example,  by  the  distribution  of  specific  enstrophy  (enstrophy  of  structure  divided 
by  the  volume  of  the  structure)  over  the  population  (Fig.  4.9). 

Specific  Enstrophy 
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Fig.  4.8:  Distribution  of  structures  across  the  height  of  the  channel  at  ReT  =  180  (Limits  of 
histograms:  0  -  360).  Left,  Tr  —  200.  Right,  Tr  =  91 
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Fig.  4.9:  Distribution  of  specific  enstrophy  of  structures  of  a  plane  channel  at  equilibrium  at 
Rer  =  180  (Limits  of  histograms:  0  -  0.3).  Left,  Tr  =  200.  Right,  Tr  =  91 

It  is  possible  to  remarked  the  much  higher  of  preferred  level  for  distribution  corresponding 
to  Tr=200  (Fig.  4.9  left)  compared  with  the  one  corresponding  to  Tr=91  (Fig.  4.9  right). 

It  can  be  concluded  that  choice  of  trigger  level  will  always  retain  a  certain  level  of  user- 
choice,  while  it  is  feasible,  and  it  might  be  useful,  to  separate  a  sub-set  of  structures  and 
study  them  in  comparison  with  the  overall  averages  of  the  whole  flow-  field. 

It  is  then  evident  that  choice  of  trigger  level  will  depend  of  the  flow  to  be  studied  and  the 
phenomenology  to  be  put  in  evidence.  In  the  case  of  attached  wall  layer,  a  strategy  would 
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be  to  take  in  account  the  different  behavior  which  are  encountered  in  the  3  layers  which  can 
be  identified  moving  away  of  the  wall 

Zone  1  Inner  (Viscous-buffer)  layer:  0  <  2+  <  30 
Zone  2  Inertial  layer:  30  <  z+  <  100 
Zone  3  Wake  region:  z+  >  100 

As  most  of  turbulent  production  is  concentrated  in  the  inertial  layer,  physics  suggests 
that  trigger  level  should  be  chosen  to  capture  most  of  it  process;  concentration  of  structures 
of  high  specific  enstrophy  (Fig.  4.9)  confirms  this  point.  Therefore,  a  general  procedure  could 
be  to  start  from  the  PDF  for  Zone  2  (Fig.  4.5)  and  choose  a  value  Tr  corresponding  to  a 
decrease  of  value  of  PDF  by  a  decade.  This  choice  is  found  to  eliminate  numerical  noise  due 
of  the  LES  algorithm  and  yet  retain  the  maximum  of  useful  information.  However,  since  the 
choice  of  this  region  retains  somewhat  arbitrary  arbitrariness,  it  was  thought  safer  to  app  y 
this  criterion  to  the  global  PDF  of  the  whole  channel,  and  this  procedure  was  scrupulously 
followed  whenever  possible,  in  the  following. 

In  present  case,  this  criterion  would  will  produce  a  value  Tr  =  91,  which  is  the  one  used 
for  the  left  sides  of  Fig.s  4.8  and  4.9. 

Structures  found  applying  this  trigger  level  have  a  content  of  enstrophy  which  is  can  be 
seen  in  Tab.  4.2  for  the  overall  channel: 


Enstrophy x 

Enstrophyy 

Enstrophyz 

Enstrophy 

47.168 

16.436 

31.859 

27.273 

Table  4.2:  Plane  channel  at  equilibrium  at  ReT  =  180  (Example).  Overall  structures  per¬ 
centage  in  enstrophy 

while,  detail  of  Zone  2,  which  contains  most  of  the  activity,  is  presented  in  Tab.  4.3. 


Enstrophyx 

Enstrophyy 

Enstrophyz 

Enstrophy 

82.882 

63.080 

64.252 

69.812 

Table  4.3:  Plane  channel  at  equilibrium  at  Rer  =  180.  (Example)  Zone  2  structures  percent¬ 
age  in  enstrophy 

showing  the  capability  this  procedure  to  isolate  the  relevant  physical  quantities  for  further 
investigation. 

To  conclude,  a  general  procedure  to  select  a  suitable  trigger  level  for  a  specific  flow  could 
be 
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1.  Select  the  region  of  interest  in  computational  field,  such  as:  inner  layer  for  attached 
wall-layer;  separate  region  or  re-attachment  region  in  wall-layers  with  detachment;  wake 
region  for  flow  around  a  body  in  free  flow;  or,  like  in  present  study,  the  overall  compu¬ 
tational  field. 

2.  Find  the  PDF  for  Q  (or  for  A2,  D ) 

3.  Fix  as  TV  an  a  value  of  Q  one  order  magnitude  lower  than  the  one  corresponding  to 
numerical  noise  of  LES 

This  process  has  been  followed  in  the  following  investigation. 

4.4  Definition  and  classification  procedure 

Once  applied  the  trigger  to  the  test  variable  (be  it  Q,  A2,  D ),  the  further  step  is  to  identify 
and  classify  the  actual  structures. 

The  final  product  of  the  detection  algorithm  is  a  "mask"  ( M(i,j ,  k)),  indicating  whether 
a  node  (M(i,  j,  k )  =  1)  satisfies  the  chosen  criterion  and  is,  potentially,  part  of  a  vortex 
structure.  All  the  grid  points  where  M{i,j,k )  =  1  are  tested  to  assess  whether  they  are 
truly  part  of  a  structure.  Basic  principle  of  acceptation  criterion  is  that  the  node  should  be 
part  of  a  connected  region  larger  than  a  single  cell,  and  the  actual  minimum  accepted  size 
is  an  octahedron  whose  vertices  are  pressure  nodes  (Fig  4.10).  All  nodes  found  within  such 
structure  or  connected  to  it  by  a  connectivity  are  tentatively  added  to  it. 


Fig.  4.10:  Smallest  Recognized  Structure 


4.4.1  Algorithm  Description 

Once  the  trigger  level  has  been  chosen,  it  is  possible  to  produce  the  mask  M(i,j,k ); 
this  mask  is  processed  looking  for  all  the  possible  nodes  that  could  be  at  the  earth  of  an 
octahedron  (Fig.  4.10).  At  this  point,  the  flag  of  interior  nodes  is  attributed  to  these  nodes 


38 


(Fig.  4.11).  The  next  logical  step  would  be  to  look  around  each  of  these  interior  nodes  to  see 
whether  the  neighboring  nodes  are  interior  as  well. 


t^v ^.t- /■> /-ti- TTndpf.prt’.ed 


Fig.  4.11:  Finding  interior  nodes 


local  surrounding  based 
onorthogonal  directions 


Fig.  4.12:  Vortex  core  identification  and  construction 

If  this  is  the  case,  and  they  are  strongly  connected,  they  are  assembled  together  giving 
birth  to  what  will  be  hereafter  called  vortex  core  (Fig.  4.12). 
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Fig.  4.13:  Skin  growth  simplified  description 


Once  the  different  vortex  cores  have  been  identified,  the  next  step  is  to  find  their  outer 
skin.  In  this  phase,  adjacent  nodes  are  checked  and  their  connectivity  tested  together  with 
neighboring  cores  in  order  to  identify  if  they  belong  to  one,  two  or  more  structures.  In  this 
phase  conflicts  could  arise,  since  the  procedure  is  typically  iterative,  meaning  only  one  node 
is  examined  at  a  time,  and  it  is  impossible  to  foreseen  whether  this  node  will  finally  belong 
to  only  one  structure  that  is  made  of  two  vortex  cores  connected  by  a  thin  membrane  or  is 
a  point  common  two  otherwise  distinct  structures. 

In  summary,  it  was  found  that  the  procedure  to  add  to  vortex  cores  a  growing  skin 
layer  is  not  enough  to  define  the  structures  in  unequivocal  way  because,  sweeping  the  field 
in  sequence,  a  node  could  be  assigned  to  different  successive  ones. 

Therefore,  it  was  found  necessary  to  iterate  over  the  flow-field,  building  virtual  list  of 
association  in  order  to  detect  all  the  conflicts  of  attribution  and  solve  them.  In  a  simplified 
way,  the  procedure  of  skin  reconstruction  can  be  resumed  with  scheme  proposed  in  Fig.  4.13. 

j  Since  the  procedure  so  outlined  is  iterative  in  a  k,j,i  fashion,  once  iterated  on  the  flow 

field,  a  list  of  possible  conflicts  is  outlined,  and  virtual  structures  grown  on  around  the  vortex 
cores  are  considered.  However,  due  to  the  iterative  procedure,  a  point  that  was  considered 
only  adjacent  to  a  structure  but  not  belonging  to  it  when  it  was  first  tested,  could,  at  the  end 
of  the  whole  scan  of  the  interested  volume,  be  indeed  found  part  of  a  structure  or  part  of  a 
filament  connecting  two  structures.  This  shortcoming  of  the  procedure  is  inherited  from  the 
procedural  language  used  to  develop  the  algorithm,  FORTRAN90,  in  the  present  case.  It  is 
not  difficult  to  imagine,  considering  the  entities  we  are  dealing  with,  that  an  object  oriented 
language  could  manage  this  phase  in  a  more  efficient  and  reasonable  way.  At  the  moment, 
the  limit  of  the  procedural  language  requires  to  reiterate  the  whole  procedure,  building  a 
library  of  virtual  structures  until  all  the  possibilities  are  examined.  In  other  terms,  a  further 
step  is  computed  to  check  whether  no  new  conflicts  are  added  to  the  previously  discovered 

f 
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in  the  first  sweep  of  the  volume  (see  Fig.  4.14). 


Fig.  4.14:  Identification  algorithm  log:  First  Fig.  4.15:  Identification  algonthm  log:  Sec- 
Iteration  ond  deration 

At  this  point  the  conflicts  are  analyzed.  Typical  outcome  of  a  conflict  analysis  would  be 
to  merge  two  or  different  structures,  while  it  is  possible  that  other  conflicts  will  not  solved 
since  they  consist  of  two  structures  briefly  touching  along  their  external  surfaces.  Once  the 
structures  database  is  updated,  the  procedure  starts  again,  looking  for  new  conflicts.  Usually, 
at  the  second  sub  step  the  number  of  conflicts  is  dramatically  reduced  (see  Fig.  4.15).  The 
procedure  continue  until  all  the  structures  are  finalized. 


Fig.  4.16:  Detection  mask  M  of  plane  channel  at  equilibrium,  at  Rer  =  180;  lower  half 
channel. 

The  result  of  present  procedure  can  be  observed  comparing  Fig.  4.16,  which  is  the  initial 
field  M(i,j.k),  with  Fig.  4.17  which  presents  the  structures  accepted  at  end  of  the  process, 
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Fig.  4.17:  Identified  structures  of  plane  channel  at  equilibrium,  at  ReT  =  180;  lower  half 
channel 


Fig.  4.18:  Noise  of  plane  channel  at  equilibrium,  at  Rer  =  180;  lower  half  channel 


while  Fig.  4.18  shows  the  nodes  accepted  from  the  detection  algorithm  and  rejected  from  the 
identification  algorithm,  or,  in  other  terms,  the  difference  between  Fig.  4.17  and  Fig.  4.16. 

For  an  acceptable  result,  the  ensemble  of  grid  cells  so  rejected  should,  ideally,  be  far 
smaller  than  the  final  retained  one;  for  the  current  case  the  corresponding  values  are  presented 
in  Tab.  4.4  and  show  that  this  requirement  is  satisfied. 

The  ratio  structures  versus  noise  could,  of  course,  be  improved  fitting  trigger  level  case  by 
case,  but  at  the  cost  the  generality  of  the  procedure;  the  present  result  shows  that  application 
of  a  method  fixed  a  priori  produces  satisfactory  filtering  of  M  and  definition  of  organized 
structures  in  flow-field,  Therefore,  it  will  be  applied  in  the  following  investigation  and  is 
recommended  for  future  applications. 

As  a  last  example,  Fig.  4.19  shows  two  rendering  of  the  same  data-field  (from  flow 
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Table  4.4:  Plane  channel  at  equilibrium  at  ReT  =  180.  (Example)  Percentage  of  Volume, 
Enstrophy,  Enstrophy*  found  in  noise  and  structures  versus  total  quantities;  entire  compu¬ 


tational  field 


Volume 

Enstrophy 

Enstrophyx 

Noise 

1.4360 

2.1259 

2.5761 

Structures 

9.8888 

27.273 

47.168 

realization  different  of  the  one  presented  in  the  previous  pages).  On  the  upper  figure  a 
particular  visualization  technique  has  been  used  instead  of  the  usual  iso-surface  procedure: 
since  to  each  node  we  consider  is  associated  a  volume,  and  in  a  staggered  grid  this  volume 
is  a  parallelepiped,  where  the  mask  M  is  equal  to  1,  there  is  a  parallelepiped  in  blue,  that  is 
transparent,  so  it  is  possible  to  see  through  them  to  see  the  identified  structures  in  red.  The 
upper  image  is  less  clear  than  the  lower  one,  where  mask  and  detected  structures  are  way 
more  simpler  to  see,  but  offers  a  glimpse  to  the  data  the  identification  algorithm  has  to  deal 
with. 

4.5  Statistical  tools 

It  was  anticipated  that  field  M,  identified  by  the  detection  algorithm,  is  treated  to  extract 
the  relevant  statistics  of  the  coherent  field.  Statistical  tools  developed  for  the  task  already 
been  applied  to  produce  the  data  presented  in  previous  section,  but  are  fully  described  and 
shown  in  this  section,  using,  as  example,  results  for  the  afore-discussed  case  of  plane  channel 
at  Rer  =  180. 

4.5.1  Geometry  and  distribution  of  structures 

Besides  overall  volume  occupied  by  the  structures  (as  given  in  Tab.  4.3),  the  distribution 
of  volume  of  individual  structures  can  be  found  (Fig.  4.20). 

Distribution  of  size  of  structures,  projected  on  the  axis  x,  y,  z,  Lx,  Ly,  Lz  can  also  be 
extracted;  at  title  of  example,  length  in  x  direction  is  presented  (Fig.  4.21). 

This  information  is  completed  with  the  position  of  the  center  of  mass  of  the  structure 
(as  already  seen  in  Fig.  4.8),  which  allows  to  localize  the  regions  of  highest  activity. 

4.5.2  Dynamics  of  structures 

Dynamics  of  structures  and  their  interaction  with  background  are  described  by 


•  Kinetic  Energy 
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Fig.  4.19:  Detected  Mask  (blue)  and  Identified  Structures  (red)  at  comparison.  Lower  half 
of  the  plane  channel  at  Rer  =  180.  View  from  top. 


KE  = 


JU^Wfdxdydz 

D 
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Fig.  4.20:  Structures  statistics:  volume  in  wall  units  Rer  (Example);  inertial  layer  (30  < 
z+  <  100).  (Limits  of  histograms:  0-7) 


Fig.  4.21:  Structures  statistics:  Stream-wise  extension  at  Rer  (Example)  in  wall  units;  inertial 
layer  (30  <  z+  <  100).  (Limits  of  histograms:  0  -  1500) 


Turbulent  Kinetic  Energy 


TKE 


=  ///> 


’dxdydz 


Vorticity, component  per  component,  here  shown  for  stream-wise  direction 

ux  =  JJJ  (V  A  v)x  dxdydz 
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•  Vorticity,  modulus 


•  Enstrophy 


■III”  A  F|  dxdydz 


=  JJJ  ^  (V  A  v)2  dxdydz 


Enstrophy,  component  per  component,  here  shown  for  stream-wise  direction 


=  JJJ  ^(V  A  v)2x  dxdydz 


where  the  Q  represents  the  whole  computational  domain;  to  all  respects,  these  are  global 
quantities. 

The  quantitative  importance  of  structures  on  overall  flow  is  expressed  as  the  ratio  of 
structure  content  versus  entire  flow,  which,  for  general  scalar  quantity  F,  is  expressed  in 
eqn.  4.12. 


IIIstrs  F dxdydz 
If  In  F dxdydz 


(4.12) 


where  STRS  represents  the  portion  of  the  domain  occupied  from  the  various  structures  and 
as  before  is  the  whole  computational  domain. 

For  example,  the  ratios  for  module  vorticity  and  enstrophy  are  represented  in  eqn.  4.13 
and  4.14. 


I II strs  v  A  V  dxdVdz 
fffn  V  A  V  dxdydz 


(4.13) 


IIIsTRS  V  A  V  dxdVdz 

Illn  ^  ^  dxdydz 


(4.14) 


It  has  to  be  remarked  that  individual  components  of  vorticity  and  enstrophy  vectors  are 
also  examined,  to  take  in  account  the  anisotropic  nature  of  turbulence  in  shear  layer. 
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Table  4.5:  Plane  channel  at  equilibrium  at  Rer  -  180.  Percentage  of  Volume,  K.E.,  T.K.E. 


Structures 

Volume 

9.8888 

K.E. 

8.0895 

T.K.E. 

14.280 

V  or  ticity  Module 

20.850 

Enstrophyx 

47.168 

EnstrophyY 

16.436 

Enstrophyz 

31.859 

Enstrophy 

27.273 

As  already  shown  in  previous  section,  percentage  of  overall  contribution  of  structures 
to  entire  flow-field  is  available.  To  complete  the  example  presented  in  previous  section, 
corresponding  data  is  presented  in  Tab.  4.5. 

While  this  information  defines  the  overall  behavior  of  structures,  the  details  of  the  this 
population  is  better  captured  introducing  the  definition  of  specific  vorticity  and  specific 
enstrophy,  as  shown  in  eqn.  4.15. 


fffsTR 

< 

t> 

2 

dxdydz 

IfJ 

rsTRdxdydz 

(4.15) 


This  definition  allow  associate  vorticity  components,  enstrophy  or  other  quantities  of  the 
structure  to  its  volume. 

This  presentation  allows  to  put  in  evidence  the  existence  of  preferential  states  for  vorticity 
and  enstrophy  associated  with  the  structures;  for  example,  specific  vorticity  and  specific 
enstrophy  of  x  axis  are  presented  in  Fig.  4.22  and  4.23 

In  previous  figures,  histograms  have  been  applied  in  order  to  approximate  numerically 
PDFs.  In  the  same  way  is  possible  to  approximate  a  joint  PDF,  which,  assuming  the  phe¬ 
nomenon  to  obey  to  the  ergodic  principle,  expresses  the  probability  for  two  variables  x,  y  of 
falling  in  a  certain  interval  centered  around  two  valuesX  and  Y\  calling  P{X,Y,dx,dy )  the 
aforementioned  probability,  we  find  the  JPDF  defined  in  eqn.  4.16. 

rX+\dx  rY+^dy 

P(X,Y,dx,dy)=  /  JPDF{xux2)dx1x2  (4.16) 

JX-\dx  JY-\dy 


47 


Cat*  B  (Tr  =  91)-  number  of  bln»:  200 


Fig.  4.22:  Structures  statistics  at  Rer  (Example):  specific  vorticity  modulus  (non- 
dimensional);  inertial  layer  (30  <  2+  <  100)  (Limits  of  histograms:  16  -  130). 


Ca*e  B  (Tr  a  91)-  number  of  bins:  200 


Fig.  4.23:  Structures  statistics  at  Rer  (Example):  specific  enstrophy  in  stream- wise(non- 
dimensional);  inertial  layer  (30  <  z+  <  100).  (Limits  of  histograms:  0  —  0.24) 


JPDFs  can  be  used  to  put  in  evidence  preferred  states  in  the  present  distribution,  in 
terms  of  structure  geometry  (Fig.  4.24),  or  vorticity  (Fig.  4.25).  For  the  latter,  it  is  useful 
to  note  that,  for  the  study  of  channel  flow,  a  proper  choice  of  the  post-processing  procedure 
will  be  able  to  extract  from  the  whole  channel  the  behavior  of  a  simple  half  channel.  This 
capability  will  be  used  following  analysis. 

A  last  quantity  of  fundamental  interest  to  understand  the  present  phenomenology  is  the 
convection  velocity  of  the  centroid  of  the  structures,  defined  in  eqn.  4.17. 
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JPDF  (%)  (Zona  2) 


Fig  4  24-  Structures  statistics  Rer  =  180  (Example):  JPDF  for  structures  length  in  x  and 
y  in  wall-units;  inertial  layer  (30  <  z+  <  100).  (Limits  of  histograms:  18  -  1400,  6  -  o50, 

200  x  200  bins) 


Spec.  Vort.  In  stroam-wfae 

Fig.  4.25:  Structures  statistics:  JPDF  for  specific  vorticity  in  x  and  y  (non-dimensional); 
inertial  layer  (30  <  <  100).  (Limits  of  histograms:  18  -  1400,  6  -  550,  200  x  200  bins) 


IIfsTRdxdydZ 


(4.17) 


which  would  allow  to  track  structures  in  time  and  assess  how  they  follow  the  average  flow. 

The  above  presented  tools  allow  to  study  in  depth  the  phenomenology  of  organized 
structures  and  will  be  applied  in  next  Chapters. 
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Ca«»  B  (Tr  a  91)-  number  of  bins:  200 


Fig.  4.26:  Structures  statistics  Rer  =  180  (Example):  Convection  velocity  (non-dimensional); 
inertial  layer  (30  <  z+  <  100).  (Limits  of  histograms:  9.5  -  19) 

4.5.3  Kinetic  energy  budget  for  structures 

To  complete  the  study  of  coherent  structures,  it  is  necessary  to  examine  their  contribu¬ 
tions  to  the  different  processes  involved  in  the  transport  of  kinetic  energy. 

Matter  of  fact,  as  it  was  pointed  out  in  paragraph  4.1,  it  is  of  fundamental  interest  to  find 
out  whether  turbulent  statistics  are  different  of  organized  field  and  un-coherent  background. 
For  this  purpose  the  budget  of  resolved  kinetic  energy  is  computed;  the  definitions  to  be 
applied  are 


Ui{  x,  t) 
ui{ x,  t) 

k 


Ui{x,t)+vSGSi{, 

x,t) 

(4.18) 

(Ui{x,t)  -  ( Ui{x 

'.*)» 

(4.19) 

1 

rUj  Ui 

(4.20) 

z 

(4.21) 

The  basics  of  transport  equation  for  kinetic  energy  in  LES  have  been  presented  in  Chapter 
2;  in  the  following,  resolved  turbulent  kinetic  energy  ( r.t.k.e .)  (which  is  indicated  as:  k )  will 
be  considered. 

For  a  statistically  stationary  flow  and  applying  the  Reynolds  decomposition,  it  is  possible 
to  identify  the  fluctuating  component  of  the  resolved  velocity  field  and  fro  it  compute  the 
instantaneous  r.t.k.e.. 

Manipulating  the  filtered  dimensional  Navier  Stokes  equations  and  averaging  in  time, 
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the  corresponding  budget  of  the  resolved  turbulent  kinetic  energy  ( r.t.k.e.b .)  is  defined 


0(*> 

dt 


( ui  uj) 


am 

dxj 


9(k  uj) 

dxj 


+ ... 


d{UiTij)  .  /jgduA 
dxj  \  ij  dxj 


(4.22) 

(4.23) 

(4.24) 

(4.25) 


Where 


P  mod  — 

t 

P  mod 


P  Uk 
p  3 

P  mod  (P  mod  ) 


(4.26) 

(4.27) 

(4.28) 


For  the  present  purpose,  a  conditional  averaging  must  be  performed  over  population  of 
coherent  structures  and  incoherent  turbulence.  This  conditional  sampling  would  provide  a 
physical  understanding  of  behavior  of  the  two  different  states  for  the  different  terms  of  the 
turbulent  kinetic  energy  budget. 

The  application  of  this  algorithm  to  the  present  case  of  plane  channel  at  equilibrium  of 
Rer  =  180  is  here  shown;  DNS  data  taken  from  AGARD  (1998)  are  given  as  reference. 

The  complete  budget  is  presented  in  Fig.  4.27  with  an  example  of  conditional  averaging: 
LESin  represents  the  conditional  averaging  performed  inside  the  coherent  structures,  while 
LESout  on  the  incoherent  turbulence. 

In  Fig.  4.29  and  4.28  a  comparison  of  the  production  and  of  the  dissipation  of  the  resolved 
turbulent  kinetic  energy  together  with  corresponding  DNS  data. 

It  can  be  remarked  that  the  behavior  and  value  for  the  different  terms  of  the  TKE 
budget  obtained  by  present  simulation  are  practically  identical  to  the  DNS  ones  except  for 
the  inner  wall  layer  of  distance  from  the  wall  smaller  than  20  wall  units.  The  peak  in  the 
turbulent  production  is  not  captured  exactly,  but,  in  the  LES  framework,  the  comparison 
can  be  assessed  satisfactory. 

Considering  the  turbulent  dissipation,  the  SGS  contribution,  which  represents  the  effects 
of  the  sub-grid  scale  motions  to  this  process,  coincides  with  the  DNS  dissipation  above  the 


51 


-->LES,.-->LESln,~->LES0ut 


Fig.  4.27:  Resolved  Turbulent  Kinetic  Energy  Budget  Rer  =  180  (Example) 


Blue  ->  LES,  Red  ->  LES. ,  Black  ->  LES.  . 

m  out 


Fig.  4.28:  Resolved  Turbulent  Kinetic  Energy  Pseudo  Dissipation  ReT  =  180  (Example) 


20  wall  units  limit.  Even  closer  to  the  wall  the  DNS  behavior  is  correctly  reproduced.  Note 
that  the  lower  production  predicted  by  LES  with  respect  to  LES  is  also  compensated  by  a 
lower  dissipation,  ensuring  a  correct  overall  budget. 

Concerning  the  conditional  averaging  presented  in  Fig.  4.27,  4.28  and  4.29,  the  results 
seem  suggest,  that  for  this  particular  simulation  flow,  the  overall  behavior  of  the  conditional 
averages  converges  to  the  classical  one  in  the  long  term,  consistent  with  Hussain  (1983) 
interpretation.  Further  tests,  performed  in  this  frame,  confirmed  this  interpretation 
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5.  SIMULATION  OF  COHERENT  STRUCTURES  IN  ATTACHED 
WALL  SHEAR  LAYER 

As  announced  in  Introduction,  the  turbulent  plane  channel  at  equilibrium  (Fig  5.1)  is 
chosen  as  test  case  for  the  study  of  coherent  structures  in  wall  shear  layers. 


x  (stream-wise) 


Fig.  5.1:  Plane  channel  geometry 

The  turbulent  flow  in  a  plane  channel  is  a  well  studied  case  and  accurate  DNS  results  at 
Reynolds  ReT  =  180, 395  and  590  are  available  (AGARD,  1998)  in  literature  for  validation. 
The  advantage  of  focusing  on  this  case  lies  with  its  geometrical  simplicity,  which  makes  it 
feasible  to  simulate  the  full  phenomenology  of  wall-bounded  flows  and  to  perform  accurate 
simulations  with  relative  limited  computer  power.  At  the  same  time,  previous  studies  (Black¬ 
burn  et  al.,  1996;  Jeong  et  al.,  1997)  have  shown  the  coherent  structures  found  in  this  flow 
to  be  fully  representative  of  general  behavior  to  be  encountered  in  this  class  of  flows  (Chacin 
and  Cantwell,  2000;  Chong  et  al.,  1998). 

Some  results  for  the  case  of  ReT  =  180  have  already  been  presented  in  Chapter  4,  as 
title  of  example,  to  demonstrate  capability  and  working  of  present  algorithm,  but  without 
detailed  description  of  this  simulation  and  its  results,  which  will  now  be  discussed  in  present 
chapter. 

However,  it  can  be  at  once  remarked  that  this  first  calculation  puts  in  evidence  the 
existence  of  organized  structures,  elongated  in  stream-wise  direction  and  characterized  by 
clockwise  and  counter-clockwise  rotation,  as  will  be  shown  in  the  following  pages. 

Aim  of  present  Chapter  is  to  evaluate  the  procedure  outlined  in  Chapter  4,  study  its 
sensibility  to  modeling  and  computational  parameters,  and,  as  well,  to  attain  a  satisfactory 
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understanding  of  physics  of  organized  turbulence  for  this  particular  flow. 

5.1  Investigation  of  organized  structures  for  plane  channel  at  ReT  = 
180  (Case  I) 

Our  analysis  starts  from  the  study  of  a  channel  flow  at  friction  Reynolds  number  equal 
to  180  (Rer  =  180).  The  LES  results  can  be  compared  with  (AGARD,  1998)..  The  basic 
parameters  for  this  first  simulation  are  resumed  in  Tab.  5.1;  in  wall  normal  direction  a  power 
stretching  (its  coefficient  being  equal  to  1.08)  has  been  employed. 

Table  5.1:  Plane  channel  at  equilibrium  at  ReT  =  180  (Case  I).  Geometrical  Parameters 


Geometry  & 

Grid 

Test  Bed 

Lx 

Ly 

Lz 

_lZL 

2 

Nx 

Ny 

Nz 

64 

96 

64 

Ax+ 

A  y+ 

A  z+ 

22 

7.85 

1.5  <  Az+  <  16 

Smagorinky  SGS  model  with  Cs  =  0.10  and  order  accurate  discretization  in  space 
an  time  are  used  for  this  analysis;  comparison  with  reference  DNS  data  shows  that  low  order 
statistical  moments  (Fig.  5.2  and  5.3)  and  budget  of  turbulent  kinetic  energy  of  reso  ved 
turbulence  (Fig.  5.4)  show  that  the  parameters  chosen  for  this  simulation  lead  to  a  satisfactory 

simulation  of  the  flow. 

Between  the  different  detection  criteria  discussed  in  Chapter  4,  the  Q  criterion  has  been 
chosen  for  detection  of  organized  structures  in  view  of  its  computational  simplicity  and  clear- 
cut  definition  of  trigger  level,  as  it  can  be  obtained  from  PDF  of  Q  (Fig.  5.5)  applying  the 
criterion  discussed  in  Chapter  4 

Application  of  this  procedure  to  present  case  yields  value  Tr  =  72,  for  trigger  level, 
expressed  in  wall  units  non-dimensionalization. 

Statistical  information  on  structures  is  obtained  performing  ensemble  average  over  results 
from  successive  instantaneous  flow-fields,  separated  by  a  1/6  flow-through  (flow-throug 
Length  of  domain/Bulk  velocity).  Experience  has  shown  that  an  ensemble  of  more  than 
10,000  structures  are  necessary  to  ensure  reliable  statistics. 

The  overall  importance  of  organised  turbulence  respect  to  total  flow  can  be  assessed  in 
Tab.  5.2,  where  percentage  of  its  contribution  to  global  quantities  of  the  turbulent  field  are 

presented. 

It  is  to  be  remarked  that  contribution  of  organized  turbulence  to  turbulent  kinetic  en¬ 
ergy  and  vorticity  is  much  higher  than  the  volume  occupied  by  the  structures  themselves, 
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Fig.  5.2:  Mean  stream-wise  velocity  and  rms  of  resolved  turbulence  of  plane  channel  at 
Rer  =  180  (Case  I) 


confirming  they  correspond  to  the  most  active  regions  of  the  flow. 

Distribution  of  enstrophy  along  the  three  axis  (Tab.  5.3)  confirms  the  existence  of  pref¬ 
erential  pattern  of  rotation  of  around  x  axis,  which  was  already  visible  in  Fig.  4.19  in  §  4. 


Fig.  5.3:  Anisotropic  part  of  the  resolved  Reynolds  stress  tensor  of  plane  channel  at  Rer  =  180 
and  total  shear  stress  (Case  I) 


These  already  interesting  results  can  be  enriched  by  a  detailed  analysis  of  the  statistics 
of  structures  population,  which  can  be  obtained  following  the  procedure  established  in  Chap¬ 
ter  4.  As  there  indicated,  the  computational  field  is  sub-divided,  for  analysis  purpose,  in  3 
layers  in  order  of  increasing  distance  z  from  the  nearer  wall. 
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Fig.  5.5:  PDF  of  Q  plane  channel  at  Rer  =  180  (Case  I) 


Table  5.2:  Plane  channel  at  equilibrium  at  Rer  =  180  (Case  I).  Volume,  turbulent  kinetic 
energy,  module  of  vorticity  and  enstrophy  within  structures  versus  corresponding  quantities 
contained  within  the  full  computational  field _ _ 


Volume 

9.661  % 

TKE 

14.075  % 

V  orticity  M  odul  us 

19.893  % 

Enstrophy 

25.058  % 

Table  5.3:  Plane  channel  at  equilibrium  at  ReT  =  180  (Case  I).  Enstrophy  within  structures 
versus  corresponding  quantities  with  the  full  computational  field 


Enstrophy  _x 

45.9  % 

Enstrophy  _y 

14.8  % 

Enstrophy  _z 

29.4  % 

in  different  regions  of  wall  layers. 

The  average  population  is  in  the  order  of  250  structures  for  flow  field.  The  distribution 
of  volumes  of  structures  is  presented  in  Fig.  5.6,  showing  values  in  the  range  1.68e+03  to 
4.8e+06  units  in  wall  non-dimensionalization,  for  an  average  value  Volmean  =  1.48e+05  wall 
units.  In  this  figure  the  usual  distinction  introduced  above  about  distance  from  the  wall  is 
used  in  order  to  empathize  a  possible  difference  in  behavior  between  the  different  zones. 

The  knowledge  of  structures  volume  is  possibly  not  interesting  in  itself,  but,  when  used 
in  conjunction  with  data  retrieved  for  the  structures  moment  of  inertia,  could  give  interesting 
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insights  on  the  structures  dimensions  and  orientations,  more  interesting  ground-work. 


Fig.  5.6:  Structures  Volume  of  plane  channel  at  ReT  =  180  with  usual  zone  distinction  (Case 
I).  (Number  of  bins  1000;  limits  of  histogram:  0.0017  •••  4.8) 


Fig.  5.7:  Structures  mass  center  of  plane  channel  at  Rer  =  180  with  usual  zone  distinction 
(Case  I).  (Number  of  bins  200;  limits  of  histogram:  5.27  •  •  •  354.57) 

The  distribution  of  structures  mass  center  across  the  height  of  the  channel  (Fig  5.7) 
shows  that  a  large  fraction  of  structures  is  located  within  the  inertial  layer  with  a  maximum 
around  top  of  the  buffer  layer,  while  structures  are  relatively  rare  within  the  wake  region 
(Zone  3).  This  behavior  is  also  shown  by  a  side  view  of  the  field  of  detected  structures 
(Fig  5.8),  where  for  clarity  only  the  structures  of  the  lower  half  of  the  channel  are  shown  (as 
will  be  for  Fig  5.10). 

Contribution  of  organized  turbulence  for  the  different  zones  is  resumed  in  Tab.  5.4. 
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Fig.  5.8:  Side  view  of  the  field  of  detected  structures  of  plane  channel  at  ReT  =  180  (Case  I) 


Table  5.4:  Plane  channel  at  equilibrium  at  ReT  =  180  (Case  I).  Volume,  turbulent  kinetic 
energy,  module  of  vorticity  and  enstrophy  within  structures  versus  corresponding  quantities 


computed  zone  by  zone 


Zone 

Zone  1 

Zone  2 

Zone  3 

Volume 

4.858  % 

21.571  % 

1.233  % 

TKE 

4.156  % 

25.825  % 

1.107  % 

Vorticity  Modulus 

5.551  % 

41.167  % 

3.236  % 

Enstrophy 

4.965  % 

66.443  % 

6.645  % 

The  collected  data  confirm  patterns  which  might  be  assumed  from  Fig.  5.7:  most  of  the 
activity  resides  in  the  inertial  layer  and  the  contribution  of  the  wake  region  is  negligible  and 
can  be,  partly,  neglected,  as  it  will  be  done  in  following  analysis. 

Advection  velocity  of  structures  mass  center  (Fig.  5.9)  indicates  they  follow  average  flow, 
without  noticeable  lag1. 

Visualization  of  whole  field  of  detected  structures  for  present  simulation  (Fig.  5.10), 
confirms  the  tendency  of  structures  to  be  elongated  in  stream-wise  direction. 

Visualization  of  individual  structures  also  puts  in  evidence  the  large  spread  of  possible 
sizes  within  this  class  (Fig.  5.11). 

However,  structures  of  much  more  complex  geometry  can  be  encountered  (and  most 
importantly,  the  algorithm  is  ‘able’  to  recognize  them),  as  shown  in  (Fig.  5.12);  if  those 
structures  truly  exist  or  the  present  trigger  level  allows  different  structures  merge, it  is  a 
point  which  needs  to  be  assessed. 

Influence  of  grid  resolution  on  the  definition  of  the  individual  structures  will  be  addressed 
in  a  following  section;  it  can,  however,  be  remarked  that  horseshoe  vortex  are  extremely  rare 

1This  comparing  the  velocities  corresponding  to  the  peaks  of  the  PDF  in  the  figure  to  an  average  velocity 
for  the  corresponding  Zone. 
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Conviction  velocity  along  stnam-wtso  (wall-units) 


Fig.  5.9:  Stream-wise  advection  of  structures  center  of  mass  of  plane  channel  at  Re 
(Case  I).  (Number  of  bins  200;  limits  of  histogram:  4.98  •  •  •  19.33) 


Fig.  5.10:  3-D  overview  of  the  field  of  detected  structures  of  plane  channel  at  Rer  =  180 
(Case  I) 


with  current  resolution,  and  a  clear  sample  was  not  possible  to  extract. 

The  statistical  analysis  of  population  of  structures  shows  as  well  the  occasional  existence 
of  relatively  long  structures,  comparable  in  their  size  to  dimension  of  the  computational 
domain.  Possibly  these  occurrences  are  linked  to  the  complex  structures  shown  in  Fig.  5.12, 
however  genuinely-looking  long  structures  can  be  identified  as  well  (Fig.  5.13). 

The  possible  complexity  of  structures  makes  it  very  difficult  to  assess  their  dimensions 
and  orientation  for  the  whole  population;  relevant  statistics  can  be  gathered  from  the  analysis 
of  projection  the  structure:  Lx,  Ly  and  Lz  on  x,  y,  z  axis  respectively  (Fig.  5.14),  which  bound 
the  cartesian  box  containing  the  actual  structure. 


Fig.  5.11:  Typical  structures  elongated  in  stream-wise  direction  of  plane  channel  at  Rer  180 
(Case  I) 

It  is  clear  from  Fig.  5.14  that  even  with  knowledge  of  Lx,  Ly,  Lz  an  important  uncertainty 
remains  concerning  their  actual  size  and  orientation  of  the  individual  structure  (see  Fig.  5.17 

below). 

Actual  distribution  of  Lx,  Ly  and  Lz  for  present  case  is  presented  in  Fig.  5.15. 

The  distribution  of  Lz  in  Fig.  5.15  needs  some  further  comment:  present  grid  is  stretched 
in  the  wall  normal  direction,  while  it  is  uniform  in  the  remaining  directions.  Non  uniformity 
of  grid  increases  the  dispersion  of  the  data  relative  to  the  size  of  the  structures  in  wall-normal 


Fig.  5.12:  Example  of  complex  structures  of  plane  channel  at  Rer  =  180  (Case  I). 


direction.  To  give  an  idea  of  better  overall  trend,  the  lines  which  bound  the  region  where  the 
data  are  collected  are  superposed. 

Observing  Fig.  5.15,  it  can  be  remarked  that  the  distributions  of  Lx,  Ly  and  Lz  are  rela¬ 
tively  similar  for  Zones  1  and  2,  \vhich  apparently  indicates  that  structures  in  the  two  layers 
could  belong  to  a  same  phenomenology,  the  most  common  values  for  Zone  2  are  indicated  in 
Tab.  5.5. 

It  has  already  been  remarked  that  the  distribution  of  Lz  appears  very  noisy,  in  Fig.  5.15, 
compared  to  the  ones  for  Lx,  Ly:  this  behavior  is  due  to  much  refined  grid  in  2  direction  and 


Fig.  5.13:  Examples  of  elongated  structures  at  Rer  =  180  (Case  I) 


Fig.  5.14:  Estimation  of  size  of  structures  from  their  projections  on  coordinate  axis 


Table  5.5:  Plane  channel  at  equilibrium  at  ReT  =  180  (Case  I).  Most  common  values  for  L 


Ly  and  Lz  in  wall  units  (Zone  S 

0 

Lx 

Ly 

Lz 

80-90 

30-40 

30-40 

Fig.  5.15:  Structures  statistics  at  Rer  =  180  (Case  I):  length  in  wall  units.  (From  top  left 
to  bottom  right:  Number  of  bins,  2000,  8000,  2000;  limits  of  histograms:  44.17  •  •  •  1413.7, 
15.7  •  •  •  463.38,  5.48  •  •  •  205.28) 


to  its  stretching;  however  the  primary  trend  can  still  be  detected,  and  it  will  be  further  put 
in  evidence  applying  the  joint  probability  function  (JPDF)  approach  introduced  in  §  4.5.2. 

In  this  case,  Lx,  Ly  and  Lz  values  are  coupled  two  by  two,  to  obtain  the  joint  probability 
function  for  the  couple  of  variables  Lx  —  Ly,  Lx  —  Lz  and  Ly  —  Lz. 

Corresponding  distributions  for  Zones  1  and  2  are  presented  in  Fig.  5.16;  following  pre¬ 
vious  remarks,  Zone  3  has  been  neglected.  Since  we  are  interested  mainly  in  finding  the  most 
probable  couple  of  values  from  a  JPDF,  this  value  is  loosely  dependent  from  the  number  of 
bins  and  the  limits  of  the  histograms  used  to  build  the  JPDF;  so,  in  order  not  to  clutter 
furthermore  the  information  provided,  the  relative  parameters  are  not  included  in  the  figure. 

It  can  be  again  remarked  that  the  main  trends  are  qualitatively  similar  for  the  first  two 
layers  (see  Tab.  5.6). 


More  probable  value  Lz  is  higher  in  Zone  2  than  one  for  Zone  1  could  be  thought  of  an 


Fig.  5.16:  Structures  statistics:  JPDF  for  structures  extension.  ReT  =  180  (Case  I).  Left 
column  Zone  1,  right  column  Zone  2,  xy- top,  xz- middle  and  yz-bottom 


indication  of  tilting  of  the  structures  away  the  from  x-y  plane  with  increasing  distance  from 
the  wall,  or  of  an  effect  of  grid  stretching. 

The  information  from  JPDF  could  be  applied  to  build  a  model  of  typical  structures: 
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Table  5.6:  Plane  channel  at  equilibrium  at  Rer  =  180  (Case  I).  Most  common  values  for  joint 
Lx-Ly ,  Lx-Lz  and  Ly-Lz  in  wall  units. _ 


Zone 

Sl 

1 

Lx  L% 

Ly  LZ 

1 

88-33 

88-20 

33-25 

2 

88-33 

88-45 

38-44 

JPDF  Lx-Ly  could  be  thought  as  indication  of  tilt  angle  in  span-wise  direction  and  JPDF 
Lx-Lz  as  an  incidence  angle  in  wall-normal  direction.  However,  this  information  might  be 
untrustworthy,  considering  Fig.  5.17,  and  a  better  approach  to  the  problem  will  be  proposed 
in  the  following. 


Fig.  5.17:  Estimation  of  size  of  structures  from  their  projections  on  coordinate  axis  could  be 
problematic 

A  better  estimation  of  size  and  orientation  of  a  structure  can  be  obtained  from  its 
principal  inertial  moments.  These  quantities  are  computed  as  an  eigenvalues  problem  posed 
by  the  inertial  tensor  of  the  structure,  where  axial  ,  centrifugal  and  tensor  of  moment  of 
inertia  are  presented  in  eqn.  5.1. 


Ix=  vV  +  z2  p  dx  d y  d z 

JSTR 


—  1 

xy  p  dx  d y  d z 

JSTR 

f  U 

Ixy 

-/ 

Axz 

L  -  \ 

'  —T 

1  Jyx 

4 

~Iyz  ' 

1 

{  -/„ 

~Izy 

h 

(5.1) 


In  first  instance,  the  problem  is  to  find  the  eigenvalues  and  eigenvectors  of  the  inertial 
tensor,  or,  rather,  the  directions  associated  to  them.  Assuming  the  structures  to  have,  in 
general,  a  preferential  direction  (and  Fig.  4.19,  5.8,  5.10  etc.  seem  to  confirm  that),  and 
considering  the  eigenvector  associated  with  the  smallest  eigenvalue  (or,  in  other  words,  the 
axis  roughly  in  line  with  the  structure),  then  it  is  possible  to  have  an  estimation  of  the 
orientation  of  the  structures  with  better  accuracy  than  the  one  yielded  by  other  methods. 
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Fig.  5.18  produces  a  comprehensive  visual  representation  of  the  alignment  in  space  for  the 
structures;  for  each  structure  the  three  directions  associated  to  the  inertial  tensor  intersect  a 
sphere  of  radius  one  in  2  points  each,  an  entrance  point  and  an  exit  point;  in  color  code,  black 
corresponds  to  the  the  direction  of  highest  interest,  or  the  main  direction  of  the  structure.  In 
present  representation,  6  symbols  correspond  to  each  structure,  but,  in  spite  of  this  apparent 
complexity,  this  was  found  to  be  the  best  visual  representation  of  the  problem  at  hand. 

To  help  to  better  visualize  the  spatial  alignment  of  the  structures,  the  sphere  has  been 
presented  under  three  points  of  view.  On  the  figure’s  left  column  there  are  the  views  relatives 
to  Zone  1,  while  the  ones  for  Zone  2  are  on  the  right  column.  Zone  3  has  been  omitted 
altogether,  and  a  particular  handling  of  the  data  (see  Chapter  4)  has  been  performed  in 
order  to  represent  in  Fig.  5.18  the  behavior  of  one  single  half-channel,  namely  the  lower  half 
channel.  It  is  understandable,  indeed,  that  for  the  whole  channel  the  main  alignment  of  the 
structures  should  be,  on  average,  in  the  stream-wise  direction. 

Processing  the  data  in  a  this  way  would  be  indeed  this  result:  no  preferential  direction 
in  the  x-y  plane,  nor  in  the  x-z.  While  this  is  unavoidable  for  the  x-y  plane  ,  considering 
only  the  data  from  the  lower  half  channel,  or  manipulating  in  a  proper  way  the  data  from 
the  whole  channel,  it  is  still  possible  to  obtain  the  preferential  alignment  for  the  structures 
in  x  —  z  plane,  which  will  be  discussed  further  on  (see.  Fig.  5.20). 

It  is  necessary  to  mention  that  to  obtain  Fig.  5.18  it  is  necessary  to  process  a  relatively 
high  number  of  structures,  which  are  obtained  lumping  together  information  from  more  than 
one  flow  realization;  for  the  left  side  (or  Zone  1)  data  from  an  average  of  30  instantaneous 
flow-fields  is  used,  while  for  the  right  (or  Zone  2)  only  6  fields  are  used.  The  reason  of  this 
difference  lies  within  the  different  number  of  structures  which  are  present  in  the  two  zones;  as 
shown  by  previous  results,  Zone  2  contains,  by  far,  the  largest  part  of  the  population  for  the 
case  under  examination.  Fig.  5.18  has  exclusively  an  illustrative  aim,  and  where  not  properly 
generated  can  be  misleading,  as  will  be  seen  in  §  5.2.  Here  it  is  sufficient  to  confirm  the  main 
stream- wise  orientation  of  the  structures. 

Extracting  in  a  quantitative  way  the  tilt  and  incidence  angle  for  the  population  from 
these  data  yields  the  results  presented  in  Fig.  5.19  and  5.20. 

Fig.  5.19  indicates  that  most  of  the  structures  in  the  different  zones  have  a  tilt  angle  of 
less  than  10  degrees,  and  predominant  tendency  is  for  the  structure  to  lie  on  the  x  axis. 

On  the  other  side,  Fig.  5.20,  shows  average  incidence  angle  of  about  7  degrees  for  Zone 
1  and  Zone  2.  However,  while  for  Zone  1  the  peak  is  clearly  defined,  for  Zone  2  the  next 
possible  angle  considered  from  the  histogram  is  of  14  degrees  and  corresponds  to  88%  of 
the  value  of  the  peak  for  7  degrees.  It  will  not  come  as  a  surprise  then  that  in  S  5.2  the 
peak  will  be  indeed  for  14  degrees.  Similar  behavior  has  already  been  remarked  in  different 
study  (Chacin  and  Cantwell,  2000). 


Structure  mean  alignment.  Zone  1  left  column,  Zone  2  right  column 
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The  inertial  tensor  £  is  useful  not  only  to  provide  an  estimation  of  the  orientation  of 
the  structures  through  its  eigenvectors,  but  through  its  eigenvalues,  it  can  be  made  to  yield 
further  information. 


i, 
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Irwrzla  Momentuums  Driven  Date 


Fig.  5.19:  Structures  statistics  ReT  =  180  (Case  I):  Tilting  angle  extracted  from  inertial 
tensor.  (Number  of  bins:  25;  limits  of  histograms:  -90  •  •  •  90) 


Inerzla  Momentuums  Driven  Oat* 


Fig.  5.20:  Structures  statistics  ReT  =  180  (Case  I):  Incidence  angle  extracted  from  inertial 
tensor.  (Number  of  bins:  25;  limits  of  histograms:  -90  •  •  •  90) 


In  fact,  it  is  possible  to  associate  to  a  structure  an  equivalent  cylinder  of  length  l  and 
radius  r,  having  the  same  principal  inertial  moments. 

Calling  Ii  the  ith  principal  inertial  moment  in  i  direction  and  assuming  Ix  <  h  <  h,  it 
is  found  that,  in  general,  I2  »  h  »  A,  so  from  an  average  of  the  two  biggest  moments  of 
inertia  and  the  volume  of  the  structure  is  possible  to  compute  the  dimensions  of  an  equivalent 
cylinder. 

According  to  Fig.  5.21  and  5.22,  the  most  probable  values  would  be  of  90  wall-units  for 
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Intrzia  Momentuums  OrK«n  0*ti:  Cyflnd«r  approach 


Fig.  5.21:  Structures  statistics  Rer  =  180  (Case  I):  Equivalent  cylinder  radius.  (Number  of 
bins:  25;  limits  of  histograms:  4.5  •••  82.4) 


Fig.  5.22:  Structures  statistics  ReT  =  180  (Case  I):  Equivalent  cylinder  length.  (Number  of 
bins:  25;  limits  of  histograms:  21  •  •  •  815) 


the  height  of  the  cylinder  and  16  wall-units  for  its  radius.  Computing  the  equivalent  volume 
7rr2/  and  comparing  it  with  volume  obtained  from  Fig.  5.6  it  can  be  seen  in  Fig.  5.23,  that 
this  model  yields  volumes  20  —  30  %  lower  than  the  original  ones. 

The  sensibility  of  the  result  to  the  choice  of  the  equivalent  geometrical  model  for  the 
structure  can  be  appreciated  comparing  the  previous  result  to  the  one  which  could  be  obtained 
assuming  the  typical  structure  to  have  the  shape  of  an  ellipsoid  of  volume  §7ra&c,  where 
a  >  b  >  c  are  the  semi-length  on  the  3  axis  of  the  ellipsoid  or  length  of  the  three  semi-axes. 
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Incrzta  Momantuums  Orhren  Oata:  Cylinder  approach 


Fig.  5.23:  Structures  statistics  Rer  =  180  (Case  I):  Equivalent  cylinder  volume  deviation 

An  advantage  of  this  second  model  is  that,  being  the  three  semi-axes  distinct,  it  is  no 
more  necessary  to  average  J2  and  J3>  decreasing  the  level  of  arbitrariness  of  the  procedure. 

Inara  la  Momantuums  Orfran  Oata:  Ellipsoid  approach 


Fig.  5.24:  Structures  statistics  Rer  =  180  (Case  I):  Equivalent  ellipsoid  minor  semi-axis 

According  to  Fig.  5.24  to  5.26,  the  most  probable  state,  for  Zone  2,  would  be  an  ellipsoid 
having  a  2  45,  b  S  17  and  c  ^  12  in  wall  units.  It  is  possible  to  compare  the  volume  of 
the  ellipsoid  with  the  volume  of  the  corresponding  structure,  being  able  finally  to  retrieve 

Fig.  5.27. 

The  differences  between  the  two  models  can  be  appreciated  from  Fig.  5.23  and  5.27. 

The  difference  is  in  order  of  20%,  which  shows  the  difficulty  to  fit  in  a  single  model  shape 
all  present  structures,  some  examples  of  which  have  previously  shown.  It  has  to  be  remarked 


Fig.  5.25:  Structures  statistics  Rer  =  180  (Case  I):  Equivalent  ellipsoid  middle  semi-axis 


Fig.  5.26:  Structures  statistics  Rer  =  180  (Case  I):  Equivalent  ellipsoid  major  semi-axis 

f 

[' 

that  the  result  can  be  influenced  by  anomalous  structures  of  very  big  size,  such  as  the  ones 
presented  in  Fig.  5.12. 

This  behavior  could  be  due  to  the  conjunct  effect  of  insufficient  resolution  of  the  grid 
and  the  choice  of  a  too  low  Tr,  which  make  impossible  for  the  present  algorithm  to  separate 
nearby  structures;  to  address  this  problem  the  effect  of  grid  resolution  will  be  discussed  in 

|  next  section. 

i 

It  can,  however,  be  concluded  that  it  is  possible  to  produce  a  model  of  typical  structure 
which  fits,  with  reasonable  accuracy,  the  present  data  and  could  be  used  as  model  of  inter¬ 
pretation  of  the  physics.  Examination  of  actual  structures  (Fig.  5.11  and  5.12)  and  the  lesser 
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Inerzla  Momentuums  Driven  Data:  Ellipsoid  approach 


Fig.  5.27:  Structures  statistics  Rer  -  180  (Case  I):  Equivalent  ellipsoid  volume  deviation 

arbitrariness  of  the  model  (see  above)  brought  the  present  authors  to  the  conclusion  that 
the  3D  elliptic  shape  is  the  most  realistic  model  for  typical  structure  and  it  was,  therefore, 
retained  for  further  study. 

Dynamic  of  organized  vorticity  is  assessed  in  Fig.  5.28,  where  specific  vorticity  and 
enstrophy  (§4.5.2)  are  presented  component  per  component.  In  this  figure  the  percentages 
are  computed  considering  the  total  number  of  structures,  while  in  Fig.  5.29,  the  percentages 
are  computed  relatively  to  the  local  number  of  structures  present  in  the  current  zone.  This 
figure  does  not  have  all  the  curves  collapsing  the  ones  over  the  others,  however  it  shows 
a  similarity  in  the  behavior  of  Zone  1  and  Zone  2,  especially  regarding  the  vorticity  and 
enstrophy  along  the  stream- wise  direction. 

Regarding  the  vorticity  components  (Fig.  5.28  left  column),  predominant  clockwise  and 
counter-clockwise  rotation  around  the  stream-wise  axis  can  be  observed:  relatively  few  struc¬ 
tures  have  no  or  low  stream-wise  specific  vorticity;  corresponding  behavior  can  be  seen  for 
wall-normal  specific  vorticity  (although  a  loss  of  modulation  is  evident  for  Zone  2)  indicating 
structures  are  rotating  around  an  axis  in  x  —  z  plane.  It  should  be  noted  that  structures  with 
a  specific  vorticity  equal  to  zero  could  be  complex  structures  as  the  ones  depicted  in  Fig.5.12, 
where  counter  rotating  structures  are  merged  together. 

Regarding  the  specific  vorticity  in  span-wise  direction,  it  is  necessary  to  keep  in  mind 
that,  in  the  present  analysis,  the  mean  vorticity  is  subtracted  from  the  instantaneous  vorticity, 
in  the  tentative  of  removing  the  influence  of  the  shear  of  the  mean  flow.  On  the  data  collected 
under  such  conditions,  gaussian  behavior  of  specific  vorticity  in  span- wise  direction  (from 
middle  left  of  Fig.  5.28  and  5.29)  apparently  puts  in  evidence  the  lack  of  preferred  rotation 
around  this  axis;  in  reality,  in  the  two  aforementioned  figures  the  average  is  performed  on  the 
whole  channel.  It  does  not  come  as  a  surprise  the  fact  that  the  peak  of  vorticity  is  at  zero, 
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Fig.  5.28:  Structure  statistics:  Specific  Vorticity  (left)  and  enstrophy  (right)  for  x  (top), 
y  (middle)  and  z  (bottom).  Rer  =  180  (Case  I).  Percentage  computed  for  total  number 
of  structures.  (Number  of  bins:  200;  limits  of  histograms  (from  top  to  bottom),  vortic¬ 
ity:  -0.39  •••  0.43,  -1.23  •••  1.18,  -0.42  •••  0.51;  enstrophy:  0  •••  0.24,  0  •••  1.66, 


Locals  Locals  Locals 


Fig.  5.29:  Structure  statistics:  Specific  Vorticity  (left)  and  enstrophy  (right)  for  x  (top),  y 
(middle)  and  z  (bottom).  Rer  =  180  (Case  I).  Percentage  computed  zone  by  zone.  (Num¬ 
ber  of  bins:  200;  limits  of  histograms  (from  top  to  bottom),  vorticity:  -0.39  •  •  •  0.43, 
-1.23  •••  1.18,  -0.42  •••  0.51;  enstrophy:  0  •••  0.24,0  •••  1.66,0  •••  0.28) 


and  it  is  hard  to  imagine  that  the  subtraction  of  mean  vorticity  could  have  a  similar  effect. 
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Indeed,  considering  only  half  a  channel,  the  specific  vorticity  presents  a  bias  consistent  with 
the  mean  shear,  see  Fig.  5.30. 


Fig.  5.30:  Structure  statistics:  Specific  span-wise  Vorticity  for  half  a  channel;  left,  percentage 
compute  respect  the  whole  channel,  right,  percentage  computed  zone  by  zone 

Regarding  the  enstrophy  components  (Fig.  5.28  right  column),  specific  enstrophy  con¬ 
firms  the  tendency  of  the  structures  in  Zone  2  to  be  more  dense  around  a  well  defined  value; 
this  pattern  is  specially  evident  for  stream-wise  component:  this  behavior  underlines  the 
well  defined  rotation  around  this  axis.  The  wide  spread  of  values  for  span- wise  and  wall- 
normal  component  in  Zone  1  puts  in  evidence  the  transition  from  longitudinal  streaks  to 
more  complex  and  three-dimensional  organization. 

The  JPDF  of  specific  vorticity  can  be  used  to  obtain  a  more  complete  understanding 
of  the  organization  of  vorticity  field;  similarly  to  what  done  for  Fig.  5.20,  it  is  possible  to 
observe  the  behavior  of  the  whole  channel  or  only  one  half. 

For  the  whole  channel,  an  averaging  is  performed  between  the  upper  half  and  the  lower 
one,  and  is  shown  in  Fig.  5.31,  while  in  Fig.  5.32  the  data  from  both  lower  and  upper  halves 
are  rearranged  in  order  to  represent  exclusively  the  behavior  of  one  half-channel,  namely  the 
lower  one. 

JPDF  of  specific  vorticity  in  x  and  y  direction  for  Zone  1  (top  left  of  Fig.  5.31  or  Fig.  5.32) 
extends  the  conclusions  drawn  from  top  left  Fig.  5.28,  showing  the  low  probability  of  existence 
of  structures  with  no  rotation  axis  in  x—y  plane  (or,  in  other  words,  not  rotating  at  all  or  only 
rotating  along  the  wall-normal  direction),  especially  for  Zone  1;  for  structures  with  rotation 
around  x  values  of  stream-wise  vorticity  are  clustered  around  ±0.15  (Fig.  5.31,  in  units  of 
present  non-dimensionalization,  namely  lov/ut2).  Corresponding  plot  for  Zone  2  (top  right  of 
Fig.  5.32  or  Fig.  5.31)  shows  a  similar  behavior,  with  a  slightly  higher  probability  of  having 
structures  not  rotating  in  the  x  -  y  plane  2. 

2  Always  keeping  in  mind  that  the  number  of  structures  belonging  to  Zone  2  is  far  greater  of  the  one  for 
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Pig.  5.31:  Structures  statistics,  whole  channel:  JPDF  for  structures  vorticity.  Left  Zone  1, 
right  Zone  2,  ary-top,  arz-middle  and  yz-bottom.  ReT  =  180  (Case  I) 

Different  behavior  can  be  noted  examining  JPDF  of  specific  vorticity  in  x  and  z  di¬ 
rections.  Considering  both  Zones  (central  drawing  on  the  left  of  Fig.  5.31),  they  have  a 
flustering  of  points  around  (0,0);  this  is  not  a  characteristics  of  the  flow  (having  structures 


Specific  Vort Icily  In  will-norm  »l  (wall— units)  Spsclflc  Vortlcity  In  wsH-normal  (will-units) 


not  rotating  at  all  or  simply  rotating  along  span-wise  direction),  rather  a  consequence  of  the 
two  lobes  distribution  of  structures  clearly  evident  in  the  center  of  Fig.  5.32:  (0,0)  connects 
the  two  region  of  higher  density  of  structures  and  consequently  inherits  a  certain  number  o 

structures  around  it. 


Remember  that  in  Fig.  5.31,  where  the  behavior  of  the  whole  channel  is  considered,  the 
four-lobe  distribution  is  due  to  the  symmetry  of  the  channel  and  rotating  and  counter-rotating 
motions  in  the  two  wall  layers. 

Considering  again  the  behavior  of  the  single  (lower)  half  channel  (middle  left  of  Fig.  5.32), 
the  joint  PDF  clearly  identifies  the  preferred  vorticity  orientation  for  the  structures,  vortical 
structures  are  usually  streamlined  entities  developing  close  to  the  wall;  due  to  vortex  dynamics 
and  interaction  with  the  other  structures,  rotation  in  stream-wise  induces  a  well  defined 
evolution  the  other  vorticity  components,  which  is  opposite  for  the  two  halves  of  the  channel. 

Pattern  of  behavior  in  Zone  2  consistent  of  this  interpretation  can  be  observed  in  central 
drawing  on  the  right  of  Fig.  5.32:  in  this  case  it  is  clear  that  most  of  the  structures  are 
provided  of  rotation  around  the  x  and  z  axis. 

Most  probable  values  of  specific  vorticity  are  respectively  ±0.12  •  •  •  0.15  and  ±0.08  •  •  •  0.1 
for  x  and  z. 

Regarding  the  JPDF  of  specific  vorticity  in  y  and  z  direction  (bottom  of  Fig.  5.31  or 
Fig.  5.32),  the  first  remark  to  do  is  the  low  probability  ,  especially  in  Zone  1,  of  structures 
having  only  rotation  along  the  stream-wise  direction,  or  no  rotation  at  all,  meaning  that  the 
structures  are  characterized  by  off-axes  rotations.  The  results  are  characterized  by  symmetry 
respect  the  wall  normal  direction,  as  it  is  to  be  expected,  and  while  in  the  averaging  in 
the  whole  channel  produces  an  additional  symmetry  along  the  other  axis  in  Fig.  5.31,  from 
Fig.  5.32  it  is  clear  that  values  for  specific  vorticity  in  y  are  higher  than  the  corresponding 
ones  for  specific  vorticity  in  z,  even  with  subtracted  the  mean  shear,  confirming  the  tendency 
of  structures  to  tilt  away  of  the  x  —  y  plane. 

5.1.1  Conclusion 

To  conclude  the  present  analysis,  it  can  be  stated  that  the  present  data  show  the  presence 
of  population  of  organized  structures,  mostly  located  within  the  inertial  layer.  Tab.  5.4 
indicates  the  concentration  of  vorticity  and  enstrophy  with  these  structures  suggesting  the 
inevitable  conclusion  that  they  must  play  a  fundamental  role  in  formation  of  turbulence. 
Analysis  indicates  the  predominant  existence  of  structures  elongated  in  stream- wise  direction. 

Most  probable  size  seems  to  be  in  the  order  of  100  wall-units  for  the  structures  main 
axis  and  20  -  40  wall-units  for  the  two  other  axis. 

The  second  value  is  in  the  same  order  of  the  wall-normal  macro-scale,  as  obtained  from 
DNS  data,  while  the  first  one  is  considerably  lower  than  the  stream-wise  macro-scale  (180 
wall  units)  extracted  from  DNS  data. 

Fig.  5.11  allows  to  conclude  that  larger  structures  are  present  but  would  be  a  small 
minority  of  the  overall  population:  structures  with  Lx  =  180  correspond  to  4%  of  the  popu- 
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lation  in  Zone  1,  and  2%  in  Zone  2.  JPDF  (Fig.  5.16)  for  Lx  -  Lz  for  Zone  2  shows  that  all 
structures  where  150  <  Lx  <  200  fall  in  region  of  joint  probability  less  than  0.1.  Estimated 
tilting  angle  for  these  structures  lies  between  10  and  20  degrees. 

Validity  and  consistence  of  present  results  and  conclusions  must  be  verified  assessing  the 
influence  of  basic  parameters  of  simulation,  namely: 


•  Grid  resolution 

•  Detection  criterion 

•  Numerical  accuracy  of  the  solver 

•  SGS  model 

This  task  will  be  addressed  in  the  next  sections. 
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5.2  Effect  of  grid  resolution 

Effect  of  size  of  LES  grid  on  simulation  of  structures  is  assessed  in  first  instance  comparing 
the  results  presented  in  §  5.1  with  the  ones  obtained  on  different  grids  as  summarized  m 

Tab.  5.7. 


Table  5.7:  Plane  channel  at  equilibrium  at  Rer  =  180.  alternative  grid  tested  (Case  II). 
Geometrical  Parameters. 


stem _ - — 

Geometry  &  Grid 

Alt.  Grid  (Case  II) 

Lx  ^ z 

h  2 _ 

Nx  Ny  Nz 

96  128  64 

Az+  A  y+  A  z+ 

14.72  5.89  1.5  <  A^+  <  16 

It  is  to  be  remarked  that  the  present  grid  has  been  refined  both  m  stream-wise  and 
span-wise  directions,  while  the  wall  normal  resolution  has  not  been  altered,  since  it  was  ob¬ 
served  that  the  grid  is  already  sufficiently  small,  compared  to  the  typical  size  of  the  expected 
structures;  all  the  test  parameters  are  identical  to  the  previous  one  (conf.  §  5.1). 

Detailed  results  of  the  LES  computations  are  presented  in  Annex  A,  while  present  dis¬ 
cussion  focus  on  differences  between  the  two  simulations. 

Applying  established  criterion  (§  4)  for  setting  trigger  level  yields  the  value  Tr  =  91. 

Global  values  for  the  two  populations  of  structures  are  compared  in  Tab.  5.8  together 
with  their  relative  deviation  (/,*,  given  h  variable  of  interest)  respect  to  the  algebraic 
averaged  value  defined  in  eqn.  5.2.  Since  only  two  sets  of  data  have  to  be  compared 

in  present  case,  the  relative  deviations  are  equal  and  opposite  and  we  represent  in  Tab.  5.8 
only  the  absolute  percentage  value. 


/. 


1  f*  fmean 

mean  /  ,  Ji  >  <  t  dev 

n 


/, 


(5.2) 


mean 


In  Tab.  5.8  the  comparison  between  the  different  quantities  has  been  made  in  the  fol¬ 
ding  fashion:  consider  the  structures  volume;  the  volume  in  average  occupied  from  t  e 
tructures  has  been  divided  by  the  total  volume  of  the  computational  field;  the  percentage 
resence  in  volume  of  the  structures  has  been  obtained  for  the  two  grids  and  the  two  results 
hen  compared.  Another  possibility  could  be  the  direct  comparison  of  the  raw  value  of  the 
uantity  between  to  two  grids.  Possibly  this  comparison  is  stricter  than  the  former  one, 
ince  the  filtering  produced  by  the  division  for  the  total  value  of  the  quantity  under  exam 
ad  example  TKE )  might  somehow  smooth  the  grid  effect;  however  comparison  of  euienan 
uantities  (see  Tab.  5.9)  shows  them  to  be  close  enough  to  neglect  this  point. 
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Table  5.8:  Plane  channel  at  equilibrium  at  Rer  =  180.  Volume,  turbulent  kinetic  energy, 
module  of  vorticity  and  enstrophy  within  structures  versus  corresponding  quantities  within 


Grid  1 

Grid  2 

Deviation 

Volume 

9.661  % 

9.888  % 

1.2 

% 

TKE 

14.075  % 

14.280  % 

0.7 

% 

Vorticity  Modulus 

19.893  % 

20.850  % 

2.3 

% 

Enstrophy 

25.058  % 

27.273  % 

4.2 

% 

Num.of  Structures 

247 

358 

18.0 

% 

Table  5.9:  Plane  channel  at  equilibrium  at  Rer  =  180.  Total,  turbulent  kinetic  energy, 
vorticity  and  enstrophy  in  present  non-dimensional  units.  Comparison  between  Case  I  and 


Grid  1 

Grid  2 

Deviation 

TKE 

115.25 

111.67 

1.6 

% 

V  orticityM  odulus 

1457.0 

1489.4 

1.1 

% 

Enstrophy 

63603.0 

64868.0 

1.0 

% 

The  deviation  is  the  order  of  1%,  indicating  that  eulerian  field  is  sufficiently  well  con¬ 
verged. 

Concerning  the  simulation  of  the  structures,  Tab.  5.8  shows  values  for  global  quantities 
to  be  very  similar  for  the  two  cases,  apart  for  a  bigger  difference  about  the  number  of  detected 
structures;  regarding  the  second,  more  resolved,  grid,  the  flow  field  seems  to  be  characterized 
by  more  numerous  and(  smaller  structures. 

This  effect  can  be  observed,  with  some  difficulties,  in  an  overall  visualization  of  the  field, 
Fig.  5.33  and  Fig.  5.34. 

With  respect  to  the  structures  detected  on  previous  grid  (see  Fig.  5.11, Fig.  5.12  or 
Fig.  5.13),  the  number  of  detected  structures  increases,  and  most  of  the  new  ones  have 
smaller  size  in  longitudinal  direction  and  are  relatively  quasi-isotropic  ones;  an  example  of 
this  type  of  structures  is  shown  in  Fig.  5.35.  This  effect  should  not  be  surprising:  it  was 
stated  stated  in  Chapter  4  that  there  is  a  minimal  dimension  for  a  structure  to  be  detected 
and  defined.  Increasing  the  resolution,  necessarily,  leads  to  the  increase  in  the  number  of  the 

small  structures  considered. 

Structures  of  topology  comparable  to  ones  found  in  the  previous  case  (Fig.  5.11)  are 
present,  as  it  should  be  expected  (see  Fig.  5.36). 
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Fig.  5.33:  3-D  overview  of  the  field  of  detected  structures  of  plane  channel  at  Rer  -  180 
(Case  II) 
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Fig.  5.34:  Lateral  overview  of  the  field  of  detected  structures  of  plane  channel  at  ReT  180 
(Case  II) 


The  statistical  importance  of  small  structures  can  be  seen  by  the  distribution  of  volume 
of  the  population  (Fig.  5.37). 

Comparison  of  Fig.  5.6  and  Fig.  5.37  offers  a  practical  example  of  the  way  the  PDF 
here  presented  should  be  considered.  Considering  Zone  2,  the  PDF  has  a  peak  of  4.68%  for 
structures  having  a  volume  of  roughly  1.3576e+04  wall  units.  Considering  the  limits  of  the 
PDF  and  the  number  of  bins,  remembering  what  explained  in  §  4.5.2,  the  bin  extension  is 
approximately  of  3840  wall  units.  Then,  since  in  Tab.  5.8  the  average  number  of  structures 
for  field  has  been  quantified  in  358  structures,  this  means  that  there  are  -^00“  ~  9  structures 
having  a  volume  oscillating  between  1.3576e+04-3840  »  9735  and  1.3576e+  +3840  «  17416. 

Now  consider  Fig.  5.38.  In  this  figure  the  limits  of  the  histograms  are  obviously  the 
same,  because  the  data  analyzed  does  not  change,  but  the  number  of  bins  has  been  halved. 
The  peak  is  of  8.38%  for  a  volume  of  1.1656e+04  wall  units.  Now  the  bin  size  is  7680,  exactly 
the  double  of  previous  case.  Now  there  are  30  structures  having  volumes  oscillating  between 
1.1656e+04  -  7680  «  3976  and  1.1656e+04  +  7680  ~  19336. 

Now  the  interval  considered  is  twice  wider  than  the  previous  one,  but  the  number  of 
structures  included  is  three  times  bigger.  It  is  easy  to  understand  how  to  properly  compare 
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Fig.  5.35:  Increased  resolution  multiply  the  presence  of  identified  structures  of  small  dimen¬ 
sions  for  plane  channel  at  Rer  =  180  (Case  II) 

different  results  between  them  appropriate  precautions  should  be  taken3. 

If,  however,  the  main  objective  of  the  analysis  is  only  to  find  the  most  probable  di¬ 
mension  of  structures  for  two  different  flow  realizations  and  then  compare  them,  most  of 
the  restrictions  briefly  stated  can  be  disregarded  and  a  simpler  and  faster  analysis  can  be 
performed. 

3This  being  the  case  for  the  results  presented  in  Annex  G,  where  the  data  handling  took  particular 
attention  in  order  to  keep  the  comparison  consistent.  See  §  5.5  for  the  meaning  of  the  Annex 


Fig.  5.36:  Typical  structures  elongated  in  stream- wise  direction  of  plane  channel  at  Rer  -  180 
(Case  II) 

The  average  volume  (in  an  algebraic  sense)  is  now  of  Volmean  =  l-06e+  ,  roughly  50% 
smaller  than  for  previous  simulation. 

This  effect  can  be  assessed  looking  to  the  distribution  of  Lx,  Ly  and  Lz  for  this  new  grid 
(Fig.  5.39  and  Fig.  5.40). 

It  must  be  remarked  that  structures  with  Lx  in  the  order  macro-scale  (namely  the  half¬ 
width  of  the  channel)  are  still  present  (Fig.  5.36):  Fig.  5.39  top  left  shows,  on  average,  around 
5  structures  close  to  this  value;  comparing  with  the  result  of  around  9  structures  for  the  first 


Fig.  5.37:  Structures  Volume  of  plane  channel  at  Rer  =  180  (Case  II)  with  usual  zone 
distinction.  (Number  of  bins  2000;  limits  of  histogram:  135.7  •  •  •  7.68e+06) 


Fig.  5.38:  Structures  Volume  of  plane  channel  at  ReT  =  180  (Case  II)  with  usual  zone 
distinction.  (Number  of  bins  1000;  limits  of  histogram:  135.7  •••  7.68e+06) 


Table  5.10:  Plane  channel  at  equilibrium  at  ReT  =  180  (Case  II).  Most  common  values  for 
in  wall  units  for  Zone  2 


I/j  j  ly  and  L 


Lx 

Ly 

Lz 

70-60 

30-40 

30-40 

grid,  the  difference  should  not  be  statistically  important,  but  the  larger  population  of  smaller 
structures  accounts  for  an  important  change  in  relative  percentages. 


88 


Fig.  5.39:  Structures  statistics  at  ReT  =  180  (Case  II):  length  in  wall  units.  (From  top  left 
to  bottom  right:  Number  of  bins,  2000,  8000,  2000;  limits  of  histograms:  14.72  •  •  •  1413.7, 
5.89  •  •  •  553.70,  1.56  •  •  •  206.72) 


The  data  outlined  above  can  be  obtained  computing  the  histograms  in  terms  of  number 
of  structures  or  can  be  retrieved  from  the  histograms  expressed  in  percentage  multiplying  for 
the  appropriate  average  number  of  structures  for  field. 

Focusing  on  zone  2,  and  repeating  the  analysis  performed  in  §  5.1,  JPDF  for  Lx  -  Ly 
and  Lx  -  Lz  (Fig.  5.40)  are  examined  to  obtain  a  better  estimation  of  typical  structure. 


Table  5.11:  Plane  channel  at  equilibrium  at  ReT  =  180  (Case  II).  Most  probable  structure 


values  for  Lx,  Ly  and  Lz 


in  wall  units  for  Zone  2  based  on  the  JPDFs 


As  shown  in  §  5.1,  this  has  to  be  compared  to  estimation  of  size  and  orientation  of  the 
structures  as  can  be  obtained  from  their  inertial  moments. 
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Structures*  Extension  m  atream-wiaa  (wall-units) 


tt met urss'  Extension  In  stream-wtea  (wall-units) 


Structures’  Extension  In  span-wise  (wall-units) 


m 


Structures’  Extension  In  stream-wls*  (wall-units) 
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Structures  Extension*  In  strssmHwta*  (wall-units) 


Structures’  Extension  In  epan-wls*  (wall-units) 


Fig.  5.40:  Structures  statistics:  JPDF  for  structures  extension.  ReT  =  180  (Case  II).  Left 
column  Zone  1,  right  column  Zone  2,  zy-top,  rz-middle  and  yz-bottom 


Inertial  moments  of  present  population  for  Zone  1  and  2  are  presented  in  Fig.  5.41.  It 
was  considered  interesting  enough  to  present,  in  Fig.  5.42  the  alignment  for  the  structures  in 
the  third  zone,  close  the  channel  centerline.  In  this  region  the  structures  are  sufficiently  far 


away  from  the  wall  for  its  influence  beginning  to  become  negligible,  starting  a  slow  return  to 
isotropy,  as  it  clear  to  the  more  random  distribution  of  the  structures. 


Fig.  5.41:  Structures  statistics,  ReT  =  180  (Case  II):  spatial  alignment,  different  pov.  A 
Structure  mean  alignment.  Zone  1  left  column,  Zone  2  right  column 


Present  results  show  patterns  very  similar  to  the  ones  remarked  for  Grid  1.  Beside  the 
graphical  representation  of  the  structures  depicted  in  Fig.  5.41  and  Fig.  5.42,  the  statistics 
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Fig.  5.42:  Structures  statistics,  ReT  =  180  (Case  II):  spatial  alignment,  different  pov.  A 
Structure  mean  alignment.  Zone  3 


gathered  on  the  structures  lead  to  Fig.  5.43  and  Fig.  5.44. 

Comparing  the  two  set  of  information,  a  possible  incongruence  appears,  that  has  not 
been  addressed  in  §  5.1  for  simplicity,  delaying  its  study  to  now. 

From  the  graphical  representation  of  the  structures  inclination,  Fig.  5.41  top  and  bottom, 
it  is  clear  how  in  Zone  1  most  of  the  structures  are  oriented  along  the  stream-wise  direction; 
in  Zone  2,  with  increasing  distance  of  the  structures  from  the  wall, their  range  of  possible 
x  —  z  angles  increases,  while  conversely  the  one  for  the  x  —  y  apparently  decreases. 

The  observation  concerning  the  x  —  z  inclination  is  confirmed  from  the  PDF  presented 
in  Fig.  5.44  while  Fig.  5.43  does  not  indicate  a  narrowing  of  the  range  of  the  tilting  angle. 

The  reason  of  this  apparent  incongruity  lies  in  what  already  said  about  Fig.  5.18.  Matter 
of  fact,  the  population  of  structures  is  different  in  number  in  Zone  2  and  Zone  1,  with  a  much 
bigger  number  of  structures  lying  in  Zone  2.  Showing  in  the  left  and  right  columns  the 
structures  corresponding  to  the  same  number  of  flow  fields  would  have  resulted  in  a  much 


Inerzte  Momentuums  Driven  Data 


— 

— 

— 

- 1 - 

k 

-  ~ 

— 

4i*  Zone  1 
^  Zone  2 
A  Zone  3 

i 

* 

• 

« 

* 

i 

i 

i 

*. . j 

i 

* 

i 

»  / 

$ 

i 

i 

i 

f 

\ 

\ 

\ 

144 

m 

U 

— 

iA- 

X-Y  Structure  angle  in  degrees 


Fig.  5.43:  Structures  statistics  ReT  =  180  (Case  II):  Tilting  angle  extracted  from  inertial 
tensor.  (Number  of  bins:  25;  limits  of  histograms:  -90  •  •  •  90) 


Inerzia  Momentuums  Driven  Data 


Fig.  5.44:  Structures  statistics  ReT  =  180  (Case  II):  Incidence  angle  extracted  from  inertial 
tensor.  (Number  of  bins:  25;  limits  of  histograms:  -90  •  •  •  90) 


higher  number  of  structures  represented  in  the  right  column  affecting  the  readability  of  the 
data.  It  was  chosen,  as  a  rule  of  thumb,  to  roughly  show  the  same  number  of  structures  on 
left  and  right  column,  so  that  the  left  column  of  Fig.  5.41  corresponds  to  roughly  five  times 
the  realizations  (namely  flow-fields)  sampled  for  its  right  column. 

This  means  that  Zone  1  has  higher  probability  to  show,  in  Fig.  5.41,  events  of  low 
probability  (or  structures  with  high  tilting  angle)  rather  than  Zone  2,  leading  to  the  apparent 
incongruence  between  Fig.  5.41  and  Fig.  5.43:  the  former  suggests  more  structures  with  high 
tilting  angle  for  the  first  zone,  the  second  seems  indicating  the  opposite. 
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To  try  to  remedy  this  possible  source  of  confusion,  a  revised-normalized  version  of 
Fig.  5.41  and  Fig.  5.43  can  be  produced. 

In  this  case,  the  first  step  is  to  present  the  PDF  relative  to  the  structure  angles  in  terms 
of  the  percentage  computed  relatively  to  the  number  of  structures  present  in  the  relative 
zone,  and  not  to  the  total. 


Inarzk  Mom*rtuum*  Driven  Data 


Fig.  5.45:  Structures  statistics  Rer  =  180  (Case  II):  Tilting  angle  extracted  from  inertial 
tensor,  local  percentages.  (Number  of  bins:  25;  limits  of  histograms:  —90  •  •  •  90) 


Fig.  5.46:  Structures  statistics  Rer  =  180  (Case  II):  Incidence  angle  extracted  from  inertial 
tensor.  (Number  of  bins:  25;  limits  of  histograms:  —90  •  •  •  90) 

In  this  way,  an  interesting  trend  can  be  remarked,  namely  the  fact  Zone  1  and  Zone  2  are 
extremely  similar  in  the  distribution  of  the  tilting  angle  (Fig.  5.45),  while  for  the  incidence 
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angle  (Fig.  5.46)  it  is  clear  that  structures  of  Zone  1  appear  to  have  the  sharpest  defined 
preferential  angle. 

In  second  place,  it  can  be  concluded  that,  the  number  of  structures  used  in  Fig.  5.41 
being  different  (the  left  column  -Zone  1-  and  the  right  one  -Zone  2-)  they  can  be  useful  for 
illustrative  reasons  only.  Indeed,  looking  back  to  §  5.1,  no  conclusion  were  drawn  from  this 
corresponding  figure,  only  from  the  relative  PDF. 

To  correct  this,  in  Fig.  5.47  the  same  number  of  structures  has  been  used,  and  then  there 
is  no  incongruence  with  Fig.  5.45. 

Regarding  the  comparison  with  the  previous  grid,  Fig.  5.43  or  Fig.  5.45  confirms  that 
there  is  no  preferred  tilt  angle.  Existence  of  an  average  incidence  angle  (Fig.  5.44  and 
Fig.  5.46)  shows  the  tendency  of  structures  to  rotate  away  from  the  x  -  y  plane:  incidence 
angle  increases  from  Zone  1  (7  degrees)  to  Zone  2  to  reach  a  most  probable  state  of  14 
degrees;  it  will  increase  further  in  Zone  3  where  it  reaches  value  21  degrees.  Again,  it  has 
to  be  remarked  that  difference  with  results  for  Grid  1  does  not  change  substantially  the 
structures  behavior. 

Following  conclusions  derived  at  the  end  of  §  5.1,  the  typical  shape  of  the  structure  is 
assumed  to  be  an  ellipsoid  of  volume  |7ro6c,  where  a  >  b  >  c  are  the  semi-length  on  the 
3  axis,  whose  values  are  presented  in  Fig.  5.48  to  Fig.  5.50.  It  is  possible  to  compare  the 
volume  of  the  ellipsoid  with  the  volume  of  the  corresponding  structure,  being  able  finally  to 
retrieve  Fig.  5.51. 

According  to  Fig.  5.48  to  5.50,  the  most  probable  state,  for  Zone  2,  would  be  an  ellipsoid 
having  a  =  42,  b  =  14,  c  =  9  wall  units. 

From  comparison  with  corresponding  results  for  Grid  1,  these  values  confirm  the  trend 
of  more  refined  grid  to  yield  a  population  with  smaller  most  probable  sizes.  It  is  to  be 
remarked  that  the  values  2  x  a,  2  x  b  are  in  acceptable  agreement  with  the  most  probable 
state  yielded  by  the  JPDF  of  Lx-Ly-  Lz  (60  -  25  -  30)  for  Zone  2  (Tab.  5.11);  obviously, 
a  perfect  agreement  cannot  be  achieved,  since  in  Tab.  5.11  the  quantities  shown  in  Fig.  5.14 
are  considered,  and  not  the  real  dimensions  of  the  structures. 

Dynamic  of  organized  vorticity  can  be  assessed  in  Fig.  5.52. 

Distribution  of  stream- wise  component  of  vorticity  (Fig.  5.52  top)  shows  a  majority  of 
structures  to  have  a  clockwise  or  counter-clockwise  rotation  around  x  axis  in  Zone  1  and  2. 
Again,  the  distribution  is  equivalent  to  one  yielded  by  simulation  on  Grid  1. 

Wall-normal  specific  vorticity  (Fig.  5.52  bottom)  shows  a  corresponding  behavior  and 
agreement,  while  gaussian-like  behavior  of  specific  vorticity  in  span-wise  direction  (middle) 
puts  in  evidence  the  lack  of  preferential  rotation  around  this  axis;  it  must  again  be  remember 
that  for  Fig.  5.52  the  average  on  the  whole  channel  has  been  used. 
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Fig.  5.47:  Structures  statistics,  i?eT  =  180  (Case  II):  spatial  alignment,  different  pov.  A 
Structure  mean  alignment.  Same  number  of  structures  (1500)  for  Zone  1  left  column,  Zone 
2  right  column 


The  JPDF  of  specific  vorticity  for  Zone  1  and  2  are  presented  in  Fig.  5.53  and  Fig.  5.54. 
These  distributions  and  most  probable  states  are  very  similar  to  ones  presented  in  §  5.1. 


Inerzia  Momentuums  Driven  Data:  Ellipsoid  approach 


Fig.  5.48:  Structures  statistics  Rer  =  180  (Case  II):  Equivalent  ellipsoid  minor  semi-axis 


Inerzia  Momentuuma  Driven  Oata:  Ellipsoid  approach 


Fig.  5.49:  Structures  statistics  Rer  =  180  (Case  II):  Equivalent  ellipsoid  middle  semi-axis 


The  only  remarkable  difference  concerns  the  x  -  y  JPDF  in  the  proximity  of  the  wall  (Zone 
1);  comparison  with  corresponding  results  on  Grid  1  shows  that  the  "dough- nut'  distribution 
where  no  points  are  found  in  (0,0)  region  (Fig.  5.31),  is  replaced  by  a  wider  dispersion  of 
points  and,  among  the  new  structures  put  in  evidence  by  the  more  refined  grid,  no-rotating 
states  are  more  probable. 

5.2.1  Conclusion 

In  conclusion,  it  can  be  estimated  that  results  shown  in  the  present  §  5.2  lead  to  fact 
that  the  global  statistics  associated  with  organized  turbulence  are  practically  identical  for 
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Fig.  5.50:  Structures  statistics  Rer  =  180  (Case  II):  Equivalent  ellipsoid  major  semi-axis 


Inaota  Mwrwntuumj  Ort*M  Oili:  Enlpjoid  wpraach 


Fig.  5.51:  Structures  statistics  Rer  =  180  (Case  II):  Equivalent  ellipsoid  volume  deviation 

Grid  1  and  2,  which  is  also  the  case  the  case  for  the  distribution  of  vorticity  within  the 
population.  The  main  difference  between  the  two  simulations  concerns  the  size  of  structures: 
grid  refining  leads  the  zones  retained  by  Q  criterion  to  become  separated  in  smaller  entities  by 
the  detection  and  classification  algorithm.  It  can  be  suspected  that  this  element  is  critically 
affected  by  the  inner  working  and  assumptions  of  present  algorithm;  therefore  this  argument 
could  only  be  set  by  direct  simulations. 

However,  the  level  of  grid  refinement  associated  with  wall-resolved  LES  is  sufficient  to 
find  quantitative  information  of  the  fields  and  its  inner  organization  and  could  be  retained 
as  guideline  for  further  investigations. 
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Fig.  5.52:  Structure  statistics:  Specific  Vorticity  (left)  and  enstrophy  (right)  for  x  (top),  y 
(middle)  and  z  (bottom).  ReT  =  180  (Case  II).  (Number  of  bins:  200;  limits  of  histograms 
(from  top  to  bottom),  vorticity:  —0.44  •  •  •  0.43,  —0.65  •  •  •  0.55,  —0.51  •  •  •  0.64;  enstrophy: 
.2 


Fig.  5.53:  Structures  statistics,  whole  channel:  JPDF  for  structures  vorticity.  Left  Zone  1, 
right  Zone  2,  xy- top,  x^-middle  and  yz-bottom.  Rer  =  180  (Case  II) 


Specific  Vortlclty  In  well-normal  On  wall-units)  Specific  Vortlelty  In  well-normal  (In  well-units)  Specific  Vortlclty  In  span-wise  (In  wall-units) 


5.3  Detection  criterion 

Assessment  of  the  sensibility  of  present  simulation  to  the  choice  of  criterion  applied  for 
the  detection  of  organized  structures  is,  evidently,  of  paramount  importance  to  evaluate  its 
results. 

This  point  is  investigated  comparing  results  presented  in  §  5.1  and  §  5.2  with  ones 
obtained  applying  criterion  D  (see  §  4);  previous  experiences  (Giammanco,  2000)  and  lit¬ 
erature  (Chacin  and  Cantwell,  2000)  have  shown  its  results  might  be  sensibly  different  of 
the  ones  obtained  applying  Q  criterion.  On  contrary,  better  agreement  was  previously  found 
between  Q  and  A2  criteria  (Giammanco,  2000),  leading  to  the  decision  to  neglect  further 
testing  on  A2. 

In  view  of  the  sensibility  of  the  results  to  grid  refinement  (§  5.2),  D  criterion  is  applied 
on  Grid  2;  all  the  other  test  parameter  being  the  same  as  stated  in  §  5.2. 

A  main  problem  encountered  applying  of  D  criterion  lies  with  the  shape  of  its  PDF 
(Fig.  5.55):  the  lack  of  evident  plateau  of  noise  makes  it  difficult  to  apply  the  procedure 
proposed  in  §  4  to  choice  the  trigger  level.  Examination  of  Fig.  5.55  brought,  tentatively,  to 
select  the  following  trigger:  Tr  =  9.6e+06  in  current  adimensionalization;  this  value  has  been 
chosen  in  order  to  have  a  similar  percentage  of  volume  occupied  from  the  structures  to  one 
obtained  in  §  5.2. 


Fig.  5.55:  PDF  of  D  plane  channel  at  ReT  =  180  (Case  III) 

Comparison  which  Tr  values  for  Q  criterion,  obtained  applying  the  same  adimension¬ 
alization,  puts  in  evidence  the  different  numerical  working  of  the  criteria  and  the  difficulty 
to  choice  a  satisfactory  Tr  for  D  (Cantwell).  Global  values  for  the  population  of  structures 
are  presented  in  Tab.  5.12  together  with  deviation  respect  to  the  corresponding  results  for  Q 
criterion  as  explained  in  §  5.2. 
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Table  5.12:  Plane  channel  at  equilibrium  at  Rer  =  180  (Case  III).  Volume,  turbulent  kinetic 
energy,  module  of  vorticity  and  enstrophy  within  structures  versus  corresponding  quantities 


Grid  2 

D 

Q 

Deviation 

Volume 

8.646  % 

9.888  % 

6.7 

% 

TKE 

15.599  % 

14.280  % 

3.6 

% 

V  orticityM  odulus 

20.335  % 

20.850  % 

1.5 

% 

Enstrophy 

29.414  % 

27.273  % 

3.8 

% 

N umber  o  fS  tructures 

190 

358 

30.0 

% 

Good  agreement  between  D  and  Q,  concerning  for  global  quantities,  can  be  remarked, 
it  should  however  be  remembered  that  the  similarity  in  percentage  in  volume  is  a  direct 
consequence  of  the  choice  of  the  trigger;  what  it  is  interesting  instead  is  the  fact  that  a 
similar  percentage  in  volume  lead  to  a  similar  percentage  in  the  other  quantities  recalled  in 
Tab.  5.12,  meaning  that  similar  fluid  events  are  captured  from  the  two  detection  criteria. 

What  stands  out  from  the  table  as  a  major  difference  lies  with  the  number  of  de¬ 
tected  structures:  the  number  of  distinct  structures  separated  by  D  criterion  is  lower, 
with  a  deviation  from  the  average  close  to  30%  (using  eqn.  5.2),  for  an  average  volume 
of  Volmean  =  1.74712e+05  wall  units  for  the  former  and  one  of  Volmean  =  1.0596e+0S  wall 
units  for  the  latter.  This  difference  could  be  due  to  the  difficulty  to  tune  the  trigger  level  for 
D  and  to  the  inherent  difference  of  the  two  procedures. 

Since  the  grid  used  is  the  same  as  one  employed  in  §  5.2,  the  eulerian  statistics  are  the 
same,  indeed  the  same  case  has  been  processed  looking  for  structures  twice,  first  for  §  5.2 
and  now;  this  means  some  more  time  steps  have  added  to  the  averaging  procedure,  but  since 
already  in  §  5.2  previous  statistics  were  already  fully  converged,  eulerian  statistics  presented 
in  Annex  A  are  also  valid  for  the  present  discussion.  , 

Shape  of  distribution  of  the  population  is  qualitatively  similar  to  one  obtained  applying 
Q  criterion  (§  5.2),  as  shown  in  Fig.  5.56  and  Fig.  5.57. 

The  distribution  of  structures  mass  center  across  the  height  of  the  channel  (Fig.  5.59) 
again  shows  the  structures  to  be  mostly  present  within  the  inertial  layer,  with  a  maximum 
around  top  of  the  buffer  layer,  while  they  are  relatively  rare  within  the  wake  region  (Zone  3). 
Zone  2  contains  in  average,  around  67%  of  the  total  number  of  structures  contained  in  this 
field  and  is,  therefore,  the  main  focus  of  the  current  investigation. 

Advection  velocity  of  centroid  of  the  structures  (Fig.  5.60)  confirms  they  follow  the 
average  flow,  without  important  lag4. 

4 Remember  what  said  in  §  5.1 
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Fig.  5.56:  3-D  overview  of  the  field  of  detected  structures  of  plane  channel  at  ReT  -  180 
(Case  III) 


Fig.  5.57:  Lateral  overview  of  the  field  of  detected  structures  of  plane  channel  at  ReT  -  180 
(Case  III) 


Fig.  5.58:  Structures  Volume  of  plane  channel  at  Rer  =  180  (Case  III)  with  usual  zone 
distinction.  (Number  of  bins:  8000;  limits  of  histogram:  711  •••  16303) 


On  the  basis  of  results  of  §  5.1  and  §  5.2,  information  on  size  and  orientation  of  structures 
is  extracted  from  JPDF  and  inertial  moments,  even  if,  for  clarity’s  sake,  we  also  present  in 
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Structure*  Mas*  Center  (z*VQ) 


Fig.  5.59:  Structures  mass  center  of  plane  channel  at  ReT  =  180  (Case  III)  with  usual  zone 
distinction.  (Number  of  bins:  200;  limits  of  histogram:  4.15  •••  354.77) 


Convection  velocity  along  stream-wlee  (wall-units) 


Fig.  5.60:  Stream-wise  advection  of  structures  center  of  mass  of  plane  channel  at  ReT  =  180 
(Case  III).  (Number  of  bins:  200;  limits  of  histogram:  4.54  •  •  •  18.86) 


Fig.  5.61  the  distribution  of  Lx,  Ly  and  Lz. 

Table  5.13:  Plane  channel  at  equilibrium  at  ReT  =  180  (Case  III).  Most  common  values  for 
Lx,  Ly  and  Lz  in  wall  units  for  Zone  2  _ _ 


Lx 

Ly 

Lz 

60-74 

23-29 

45-47 

Inertial  moments  of  present  population 


for  Zones  1  and  2  are  presented  in  Fig.  5.63. 


It 
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Fig.  5.61:  Structures  statistics  at  Rer  =  180  (Case  III):  length  in  wall  units.  (From  top  left 
to  bottom  right:  number  of  bins:  2000,  8000,  2000;  limits  of  histograms:  29.45  •  •  •  1413, 
11.78  •  •  •  753,  4  •  •  •  190) 


was  considered  interesting  enough  to  present,  in  Fig.  5.64,  the  corresponding  results  for  Zone 
3,  close  the  channel  centerline.  In  this  region  the  structures  are  sufficiently  far  away  from  the 
wall  for  this  influence  beginning  to  become  negligible,  so  starting  a  slow  return  to  isotropy, 
as  it  is  clear  to  the  more  random  distribution  of  the  structures. 

Beside  the  graphical  representation  of  the  structures  depicted  in  Fig.  5.63  and  Fig.  5.64, 
the  statistics  gathered  on  the  structures  allow  to  define  the  distribution  of  tilt  angle  (Fig.  5.65) 
and  of  incidence  angle  (Fig.  5.66)  for  present  population. 

Fig.  5.65  confirms  that  there  is  no  preferred  tilt  angle;  variance  of  distribution  increases 
moving  from  Zone  1  to  Zone  2,  while  the  sample  data  for  Zone  3  is  too  small  to  allow 
firm  conclusions.  Existence  of  an  average  incidence  angle  (Fig.  5.66  )  shows  the  tendency 
of  structures  to  rotate  away  from  the  x  —  y  plane:  incidence  angle  slightly  increases  from 
—  7—14  degrees  in  Zone  1  to  ^  14  —  21  degrees  for  Zone  2;  agreement  with  corresponding 
results  for  Q  criterion  (§  5.2)  is  satisfactory. 
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Structures’  Extension  In  span-wise  (wall-units) 


Fig.  5.62:  Structures  statistics:  JPDF  for  structures  extension.  Rer  =  180  (Case  III).  Left 
column  Zone  1,  right  column  Zone  2,  xy- top,  xz- middle  and  yz-bottom 


The  typical  shape  of  the  structure  is,  again,  assumed  to  be  an  ellipsoid  of  volume  §7ra&c, 
where  a  >  b  >  c  are  the  semi-length  on  the  3  axis,  whose  values  are  presented  in  Fig.  5.67 
to  Fig.  5.69.  It  is  possible  to  compare  the  volume  of  the  ellipsoid  with  the  volume  of  the 
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Fig.  5.63:  Structures  statistics,  Rer  =  180  (Case  III):  spatial  alignment,  different  pov.  A 
Structure  mean  alignment.  Zone  1  left  column,  Zone  2  right  column 


corresponding  structure,  being  able  finally  to  retrieve  Fig.  5.70. 

According  to  Fig.  5.67  to  5.69,  the  most  probable  state  for  Zone  2  would  be  an  ellipsoid 
having  sizes  of:  a  =  35  -  42,  b  =  13,  c  =  9  wall  units. 

i! 

i 


Fig.  5.64:  Structures  statistics,  Rer  =  180  (Case  III):  spatial  alignment,  different  pov.  A 
Structure  mean  alignment.  Zone  3 


Values  ,  which  are,  again,  comparable  to  the  ones  yielded  by  Q  criterion,  even  if  present 
a  is  slightly  shorter  than  the  one  found  in  previous  case.  It  must  be  remarked  that  the  values 
2  x  a,  2  x  6  are  close  to  the  most  probable  state  yielded  by  the  JPDF  Lx  -  Lz  for  Zone  2: 
Lx-Lz  =  60-30  (Fig.  5.62);  this  result  confirms  the  capability  of  present  approach  to  identify 
the  dimensions  of  the  representative  structure,  in  the  limits  of  the  current  implementation  of 
detection  and  identification  algorithms. 

At  this  point  in  the  analysis,  an  apparent  contradiction  is  evident:  information  of 
Tab.  5.12  and  Fig.  5.58  suggest  that  typical  structure  defined  by  D  criterion  should  be 
much  bigger  than  the  ones  defined  by  Q  criterion  (Tab.  5.2);  on  the  contrary,  examination 
of  corresponding  JPDF  (Fig.  5.62)  and  data  from  inertial  moments  (Fig.  5.67  to  5.69)  yield 
a  typical  structure  much  smaller,  notably  equivalent  to  the  one  extracted  by  application  of 
Q  criterion  (§  5.2). 

Reason  of  this  result  lies  with  the  nature  of  the  detection  criterion  and  the  structures  so 
defined:  it  has  already  been  remarked,  in  §  5.1,  the  presence  of  abnormal  structures  of  very 
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Imrzta  Moment u urns  Driven  Data 


X-Y  Structure  angle  In  degreee 


Fig.  5.65:  Structures  statistics  ReT  =  180  (Case  III):  Tilting  angle  extracted  from  inertial 
tensor.  (Number  of  bins:  25;  limits  of  histograms:  —90  •••  90) 


Inerzia  Momentuunra  Driven  Data 


X-Z  Structure  angle  In  degreee 


Fig.  5.66:  Structures  statistics  ReT  =  180  (Case  III):  Incidence  angle  extracted  from  inertial 
tensor.  (Number  of  bins:  25;  limits  of  histograms:  —90  •••  90) 


large  sizes  (Fig.  5.12)  which  could  skew  the  interpretation  of  estimation  of  average  volume 
of  population.  Examination  of  present  population  shows,  at  same  level  of  grid  refining,  this 
behavior  is  much  more  common  among  structures  defined  by  D  criterion  than  among  ones 
yielded  by  Q  criterion. 

Examples  of  this  effect  can  be  seen  in  Fig.  5.71. 

However,  the  most  common  state  remains  made  of  structures,  of  different  sizes  (Fig.  5.72), 
elongated  in  x  direction,  similar  to  the  ones  presented  in  §  5.1  and  §  5.2. 
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Inerzia  Momantuuma  Drivan  Data:  ellipsoid  approach 


Fig.  5.67:  Structures  statistics  Rer  =  180  (Case  III):  Equivalent  ellipsoid  minor  semi-axis 

Inarzia  Momantuums  Ortvao  0*a:  GlHpaotd  approach 


Fig.  5.68:  Structures  statistics  ReT  =  180  (Case  III):  Equivalent  ellipsoid  middle  semi-axis 

It  can  be  concluded  that  the  presence  of  structures  of  abnormal  sizes  is  a  common 
problem,  related  not  only  to  the  nature  of  the  algorithm,  but,  mostly,  to  the  result  of  the 
detection  criterion  and  the  refinement  of  the  grid.  In  other  words,  the  algorithm  has  been 
thought  in  order  to  classify  connected  region  of  space  as  structures,  and,  to  be  applicable  to 
LES  resolution  (against  DNS  resolution),  it  has  been  designed  in  order  not  to  be  extremely 
strict,  leaving  open  the  possibility  of  connectivity  in  off-cardinal  directions.  Nevertheless, 
from  Fig.  5.72,  it  is  clear  that  the  algorithm  identifies  the  structures  in  a  logical  way,  in  all 
the  cases  the  present  mask  allows  it  to  do  so. 

Therefore,  definition  and  analysis  of  a  typical  structures  must  be  performed  considenng 
the  most  probable  state;  evaluating  a  ’ Volmean  of  structures  by  dividing  the  detected  volume 
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Fig.  5.69:  Structures  statistics  Rer  =  180  (Case  III):  Equivalent  ellipsoid  major  semi-axis 


Irwrzta  Momertuum*  Owen  Oats:  ElflpaokJ  appmtdi 


Fig.  5.70:  Structures  statistics  Rer  =  180  (Case  III):  Equivalent  ellipsoid  volume  deviation 

by  the  number  of  structures  found  is  proved  a  far  too  crude  approximation.  Taking  the 
data  of  Fig.  5.58  and  considering  Zone  2,  the  peak  of  the  PDF  corresponds,  roughly,  to 
a  Vol  =  1.1919e+04  in  wall-units;  comparing  this  value  to  the  most  probable  one  yielded 
applying  the  ellipsoid  model  (in  Fig.  5.73)  the  corresponding  value  for  of  Zone  2  is,  roughly, 
1.0552e+04  in  wall  units,  with  a  reasonable  agreement  (relative  error  6%)  between  the  two 
estimations. 

Furthermore,  these  values  are  very  close  to  the  ones  found  applying  Q  criterion  over  the 
same  grid;  matter  of  fact,  Vol  =  1.4739e+04  wall  units  is  found  from  Fig.  5.37,  while  from 
an  equivalent  set  of  data  for  §  5.2  (not  presented  in  present  text)  a  value  of  1.3432e+04  wall 
units  is  recovered. 
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Fig.  5.71:  Example  of  complex  structures  of  plane  channel  at  Rer  =  180  (Case  III) 


It  can  be  concluded  that  reliable  estimation  of  size  of  typical  structure  can  be  obtained 
from  inertial  moments  of  the  population,  and,  within  these  terms,  the  two  models  show  an 
acceptable  agreement. 

Dynamic  of  organized  vorticity  can  be  assessed  in  Fig.  5.74,  where  specific  vorticity  is 
presented  component  per  component. 

Distribution  of  stream- wise  component  (Fig.  5.74  top  left)  shows  a  majority  of  structures 
having  clockwise  or  counter-clockwise  rotation  around  x  axis  in  Zone  1  and  2.  A  major 
difference  versus  results  of  Q  (§  5.2),  lies  with  the  existence,  in  Zone  2,  of  an  important 
number  of  structures  having  negligible  rotation  around  the  x  axis.  This  results  underlines 
the  different  working  of  D  and  Q  criteria  and,  in  consequence,  the  different  definition  of 

structure. 

Negligible  importance  of  Zone  3  is,  again,  put  in  evidence  by  the  results.  Wall-normal 
specific  vorticity  (same  figure  bottom  left)  shows  a  corresponding  behavior,  while  Gaussian 
behavior  of  specific  vorticity  in  span- wise  direction  (middle  left)  puts  in  evidence  the  lack  of 
preferential  rotation  around  this  axis  (again,  it  must  be  considered  that  the  average  as  been 
performed  on  the  whole  channel). 

The  JPDF  of  specific  vorticity  for  Zone  1  and  2  are  presented  in  Fig.  5.75. 


113 


1400^>o 


Fig.  5.72:  Typical  structures  elongated  in  stream-wise  direction  of  plane  channel  at  Re 
(Case  III) 


Inarzia  Momtituoma  Driven  Data:  ElHpaokJ  approach 


Structrues’  volume  In  millions  of  wall-units 


Fig.  5.73:  Equivalent  ellipsoidal  Structures  Volume  of  plane  channel  at  Rer  =  180  (Case  III) 
with  usual  zone  distinction.  (Number  of  bins:  8000;  limits  of  histograms:  627  •  •  •  3.71e+  ) 


JPDF  of  specific  vorticity  in  x-y  plane  for  Zone  1  is  similar  to  the  previous  cases  (§  5.1 
and  5.2),  while  for  Zone  2,  as  a  confirmation  of  what  seen  in  Fig.  5.74,  a  certain  number 


Fig.  5.74:  Structures  statistics:  Specific  Vorticity  (left)  and  enstrophy  (right)  for  x  (top),  y 
(middle)  and  2  (bottom).  Rer  =  180  (Case  III).  (Number  of  bins:  200;  limits  of  histograms 
(from  top  to  bottom),  vorticity:  -0.26  •  •  •  0.27,  -0.56  •  •  •  0.58,  -0.44  •  •  •  0.41;  enstrophy: 
0  0.08,  0  •••  0.38,  0  •••  0.21) 
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Considering  JPDF  of  specific  vorticity  in  the  x  -  z  plane  for  Zones  1  and  2,  it  can  be 
again  remarked  that  distribution  and  value  for  most  probable  cases  are  similar  to  ones  found 
applying  Q  criterion,  with  the  exclusion  of  the  already  mentioned  structures  clustered  abound 
(0,0);  for  the  case  of  rotating  structures,  range  of  possible  values  is  similar  to  the  ones  yielded 
by  Q  criterion,  as  for  specific  vorticity  in  the  y  —  z  plane. 

In  conclusion,  application  of  Q  and  D  criteria  underlines  the  difference  between  the  two 
approaches,  already  discussed  in  §  4  and  previous  analysis  (Rambaud,  2000)  and,  notably,  the 
detection,  by  D,  of  a  larger  population  of  no-rotating  structures.  Present  results,  however, 
do  not  allow  to  definitively  recommend  one  criterion  over  another. 

In  the  opinion  of  present  authors,  one  main  drawback  of  D  criterion  lies  with  the  uncer¬ 
tainty  concerning  the  selection  of  proper  trigger  level,  while  consistent  results  for  Q  criterion 
can  always  be  obtained  applying  the  procedure  proposed  in  §  4,  as  shown  in  present  Chapter 
and  as  it  will  be  confirmed  in  chapter  6. 

Furthermore,  at  same  grid  size  and  similar  percentage  in  volume  occupied  by  the  struc¬ 
tures,  D  criterion  seems  more  likely  to  define  abnormal  structures  of  very  big  sizes,  which 
are  consequence  of  interaction  of  the  mathematics  of  the  criterion  and  present  algorithm,  at 
risk  of  skewing  the  conclusions  of  the  analysis. 

However,  it  is  has  to  be  stressed  that  analysis  based  upon  inertial  moments  of  the 
distribution  yields,  however,  similar  sizes  for  typical  structure,  while  discrepancy  of  global 
statistics  is  in  order  of  5%  and  JPDF  results  show  satisfactory  agreement  concerning  values 
and  main  trends. 

This  level  of  agreement  being  acceptable  for  present  investigation,  lower  CPU  require¬ 
ment  and  more  reliable  choice  of  trigger  level  lead  to  retain  Q  criterion  for  the  remaining  of 
the  investigation,  while  it  is  confirmed  that  inertial  moments  approach  represents  the  most 
reliable  method  to  define  size  and  orientation  of  typical  structures  for  this  flow  topology. 


5.4  Effect  of  precision  of  LES  numerical  algorithm 

As  discussed  in  §  3,  numerical  precision  of  LES  can  be  increased,  and  quality  of  results 
improved,  applying  a  formal  4t/l  order  accurate  discretization  in  space;  importance  of  this 
choice  on  the  simulation  of  coherent  structures  will  now  tested.  As  the  present  purpose  is 
only  to  evaluate  the  sensibility  of  the  simulation,  only  grid  1  (see  §  5.1)  is  tested5,  keeping  in 
mind  that  previous  conclusions  on  effect  of  grid  refining  (§  5.2)  will  also  apply  to  this  case. 
The  simulation  is  performed  applying  the  value  is  Cs  =  0.1  for  Smagorinsky  coefficient,  the 
same  used  in  §  5.1  which  was  then  found  to  yield  the  best  agreement  with  DNS  data.  In 
Annex  C,  the  eulerian  statistics  for  this  same  case  computed  for  Cs  =  0.065  are  presented 
to  allow  to  appreciate  the  influence  of  the  Smagorinsky  coefficient  on  the  low  and  high  order 
eulerian  statistics. 

Q  criterion  is  applied  and  a  trigger  level  Tr  =  88  is  found  for  this  case,  showing  an 
acceptable  agreement  with  value  found  for  test  on  2nd  order  discretization  (§  5.1). 

Deviation  for  the  global  quantities  for  the  entire  flow  field  is  shown  Tab.  5.14. 


Table  5.14:  Plane  channel  at  equilibrium  at  ReT  =  180  (Case  IV).  Total,  turbulent  kinetic 
energy,  vorticity  and  enstrophy  in  present  non-dimensional  units.  Comparison  between  2 


cretizanon _ 

Grid  1 

2nd 

4th 

Deviation 

TKE 

115.25  % 

105.93  % 

4.2  % 

V  orticity  Modulus 

1457.0  % 

1471.6  % 

0.5  % 

Enstrophy 

63603.0  % 

67477.  % 

1.0  % 

Global  statistics  of  organized  turbulence  for  the  two  simulations  are  compared  in  Tab.  5.15. 


Table  5.15:  Plane  channel  at  equilibrium  at  Rer  =  180  (Case  IV).  Volume,  turbulent  kinetic 
energy,  module  of  vorticity  and  enstrophy  within  structures  versus  corresponding  quantities 


Grid  1 

2nd 

4th 

Deviation 

Volume 

9.661  % 

7.582  % 

12.1 

% 

TKE 

14.075  % 

11.472  % 

m 

Vorticity  Modulus 

19.893  % 

16.248  % 

10.1 

% 

Enstrophy 

25.058  % 

20.000  % 

11.2 

% 

Number  of  Structures 

247 

210 

8.1 

% 

Tab.  5.14  puts  in  evidence  a  good  agreement  of  the  two  simulations. 

5Details  of  present  LES  can  be  seen  in  Annex  B  where  improved  agreement  with  DNS  data  can  be  noted. 
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On  the  contrary,  Tab.  5.15  shows  a  10%  deviation  for  the  quantities  associated  with 
the  structures,  as  it  could  have  expected,  keeping  in  mind  that  all  the  detection  criteria  are 
based  on  calculations  of  gradient  of  velocity  gradients  and,  therefore,  potentially  sensible  to 
the  effect  of  numerical  accuracy. 

However  the  main  characteristics  of  the  population  are  similar  to  the  ones  remarked 
in  §  5.1,  as  shown  in  Annex  B:  in  particular  JPDF  for  Lx,  Ly  and  Lz  show  values  and 
distribution  similar  to  the  ones  yielded  by  2nd  order  simulation  (§  5.1).  Similar  agreement  is 
also  found  for  JPDF  of  components  of  specific  vorticity. 

It  could  be  concluded  that  numerical  accuracy  of  LES  is  not  of  primary  importance, 
compared  with  the  effect  of  grid  refinement  and  detection  criterion  and  the  extra  cost  of 
this  refined  calculation  would  be  justified  only  when  very  precise  estimation  of  quantities 
associates  with  the  structures  was  necessary  for  design  purpose. 
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5.5  Effect  of  choice  of  SGS  model 

As  discussed  in  §  2,  SGS  modeling  is  the  corner  stone  of  LES,  and  its  effects  on  simulation 
of  the  structures  must,  therefore,  be  assessed6.  In  order  to  investigate  this  point,  a  comparison 
is  performed  on  grid  1  (§  5.1),  where  SGS  model  can  be  expected  to  have  a  detectable  effect 
on  the  resolved  field;  an  advantage  of  being  able  to  evaluate  the  TKE  budget  (see  §  2  and  §  4) 
lies  with  the  possibility  of  evaluating,  in  quantitative  way,  the  effect  of  the  SGS  model  in  the 
simulation:  knowing  that  in  any  case  the  flow  field  is  deeply  influenced  from  the  choice  of  the 
model  (this  effect  being  visible  in  the  resolved  terms  of  the  TKE  budget),  the  SGS  dissipation 
( eksos  in  Fig.  5.4  or  Fig.  A.3)  is  a  good  indicator  of  the  influence  of  the  SGS  model.  For  the 
simulation  considered  (§  5.1),  the  Smagorinsky  model  apparently  plays  a  small  role  (Fig.  5.4), 
the  dissipation  that  it  provides  being  quite  small  compared  to  the  viscous  dissipation,  which 
is  the  case  for  a  wall  resolved  LES  computation7,  however,  it  will  be  shown  how  the  change 
of  the  SGS  model  alters  the  statistics. 

Apart  the  already  applied  Smagorinsky  model  (lSMA’),  filtered  structure  function  (‘FSF’) 
model  (Compte  et  al.,  1994)  and  mixed  (‘SS’)  model  (self-similar  model  S3+Smagorinsky  dis¬ 
sipation,  Louda  (2001))  are  tested  together  with  un-resolved  DNS  (‘NO-SGS’)  approach. 

Results  for  the  resolved  field  for  these  simulations  are  presented  in  Annex  D  for  FSF, 
E  for  NO-SGS  and  F  for  SS.  First  of  all,  it  must  be  remarked  that,  since  the  self  similar 
model  is  not  an  eddy  viscosity  model,  a  different  version  of  the  TKE  budget  has  to  be  used; 
however,  even  if  developed,  it  has  not  been  sufficiently  tested  and  so  it  is  not  presented  in 
Annex  F. 

It  can  be  seen  that  all  the  SGS  models  are  able  provide  an  acceptable  reproduction  of 
DNS  data,  the  biggest  scattering  can  be  found  in  the  statistics  in  the  cross-flow  direction, 
starting  from  low  order  statistics  like  the  velocity  rms.  Some  interesting  remarks  can  be 
extracted  from  the  comparison  of  the  TKE  presented  in  §  5.1  and  Annex  E  (NO-SGS):  there 
is  a  over-prediction  of  rms  associated  to  span-wise  and  wall  normal  direction,  as  well  of 
turbulent  kinetic  energy  production  and  in  turn  dissipation:  the  eddy  viscosity  model  is  well 
suited  for  this  kind  of  flow.  In  Fig.  E.3  it  is  possible  to  observe  how  the  unresolved  LES  is 
still  very  close  to  the  DNS  behavior,  at  least  in  term  of  shape  for  the  different  TKE  terms. 

Now,  for  an  interesting  remark  about  the  difference  between  two  SGS  models,  consider 
first  Fig.  E.2  and  5.2:  apparently  for  the  FSF  model,  fluctuations  are  more  damped,  in  partic¬ 
ular  in  span-wise  and  wall  normal  directions,  while  the  turbulent  shear  stress  has,  apparently, 
a  better  behavior;  apparently  this  SGS  model  is  more  dissipative  than  the  Smagorinsky,  and 
this  is  clearly  visible  in  Fig.  E.3  and  5.4,  where  the  SGS  dissipation  ( ekscs )  is  now  almost 

6 Besides  the  choice  of  model,  the  tuning  of  the  empirical  constants  involved  should  furthermore  taken  into 
account,  see  §  5.4  and  difference  between  Annex  B  and  C 

t'j’his  observation  could  provide  a  nice  and  obvious  way  of  classifying  LBS  computations,  according  from 
the  relative  magnitude  of  the  viscous  dissipation  vs  SGS  dissipation 
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comparable  to  the  viscous  one  (ek).  On  the  other  side,  the  SGS  contribution  of  Smagor insky 
to  the  total  stress  is  negligible  (Fig.  5.3),  while  FSF  has  a  decisive  influence  on  the  shear 

stress  (Fig.  E.l). 

Like  stated  above,  all  the  SGS  models  are  able  provide  an  acceptable  reproduction  of 
DNS  data,  but  an  important  scattering  exists  for  the  results.  Total  values  for  global  quantities 

are  presented  in  Tab.  5.16. 


Table  5.16:  Plane  channel  at  equilibrium  at  Rer  =  180  (case  V).  Global  quantities:  turbu¬ 
lent  kinetic  energy,  vorticity  and  enstrophy  in  present  non-dimensional  units.  Comparison 


TOTAL 

NO-SGS 

SMA 

FSF 

SS 

TKE 

119.79 

115.25 

110.64 

95.16 

Vorticity  Modulus 

1558.0 

1457.0 

1286.9 

1279.8 

Enstrophy 

73131.0 

63603.0 

48161.0 

48060.0 

DEVIATION 

NO-SGS 

SMA 

r“  FSF 

SS 

f  TKE 

8.7% 

4.6  % 

0.4  % 

13.7  % 

Vorticity  Modulus 

11.7  % 

4.5  % 

7.8  % 

8.3  % 

Enstrophy 

25.6  % 

9.2  % 

17.3  % 

17.5  % 

A  typical  deviation  in  the  order  of  10%  for  turbulent  kinetic  energy  and  vorticity  is  clear 
from  Tab.  5.16. 

In  comparison,  percentages  of  contribution  of  organized  field  show  much  smaller  devia¬ 
tion,  as  shown  by  Tab.  5.17. 

Detailed  comparison  is  presented  in  Annex  G  and  shows  a  similar  agreement  in  what 
concerns  the  statistical  behavior  of  the  population.  t 

It  can  be  concluded  at  the  current  stage,  that  the  search  algorithm  is  always  able  to 
isolate  the  active  part  of  this  field  and  the  choice  of  SGS  model  must  be  based  on  quality 
of  reproduction  of  global  resolved  turbulence  or/and  the  cost /effectiveness  of  the  actual 

calculation. 
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Table  5.17:  Plane  channel  at  equilibrium  at  Rer  =  180.  Volume,  turbulent  kinetic  energy, 
vorticity  and  enstrophy  in  present  non-dimensional  units.  Comparison  between  NO-SGS, 


SMA,  FSF  and  SS  _ 

NO-SGS 

VOLUME  8.066  % 

TKE _ 12.357  % 

Vorticity  Modulus  17.666  % 

Enstrophy  23.32  % 

Num.  Structures  280 

DEVIATION 
Volume 
TKE 

Vorticity  Modulus 
Enstrophy 


SMA 


9.661  % 


14.075  % 


19.893  % 


25.06  % 


247 


FSF 


9.350  % 


12.651  % 


18.32  % 


21.217  % 


210 


SS 

8.655  % 
13.189  % 
18.929  % 
23.86  % 
231 


NO-SGS 

SMA 

FSF 

SS 

9.7 

% 

8.1 

% 

4.7 

% 

3.1 

% 

5.4 

% 

7.7 

% 

3.2 

% 

3.8 

% 

5.5 

% 

6.4 

% 

2.0 

% 

1.2 

% 

16.0 

% 

2.0 

% 

13.0 

% 

4.5 

% 
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5.6  Conclusions 

Present  results  have  demonstrated  capability  of  LES  to  reproduce  the  behavior  of  or¬ 
ganized  turbulence  and  ability  of  present  algorithms  to  detect  and  classify  the  recognized 

entities  embedded  in  the  flow-field. 

Structures  population  statistics  and  their  relative  importance  with  respect  to  the  total 
turbulent  field  can  be  put  into  evidence. 

Global  statistics  of  the  field  of  coherent  structures  were  found  to  be  relatively  insensible 
to  different  simulation  parameters  tested,  with  the  exception  of  what  concerns  shape  and  size 
of  individual  structures,  which  is  critically  dependent  on  grid  size  and  detection  criterion. 

This  last  result  had  to  be  expected,  in  view  of  the  fact  this  element  is,  necessarily, 
strongly  dependent  of  the  numerical  detail  of  the  implementation. 

As  this  point  is  not  of  paramount  importance  for  the  understanding  of  the  physical  be¬ 
havior  of  the  overall  phenomenon  it  is  decided,  and  is  recommended  for  future  investigations, 
on  ground  of  cost-effectiveness,  to  prosecute  the  study  applying. 


•  Smagorinsky  model  as  SGS  model 

•  2nd  order  accuracy  of  space  discretization 

•  Q  criterion  of  detection  of  structures 
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6.  INVESTIGATION  OF  COHERENT  STRUCTURES  AT  DIF¬ 
FERENT  REYNOLDS  NUMBERS 

In  previous  chapters,  phenomenology  of  organized  turbulence  was  studied  for  the  case 
of  channel  flow  at  Rer  =  180,  in  order  to  limit  the  computational  effort  involved  with  the 
assessment  of  performances  of  present  algorithm  and  its  sensibility  to  different  parameters  of 
simulation. 

Regardless  from  the  data  until  now  gathered,  results  and  conclusions  reached  cannot  be 
considered  fully  representative  of  behavior  and  characteristics  of  coherent  structures  in  fully 
developed  turbulence,  since  the  Reynolds  number  adopted  was  low  enough  for  flow  to  be 
affected  by  viscous  effects. 

In  order  to  remove  this  deficiency  and  complete  the  present  study,  flow  at  Rer  =  395 
and  590,  where  turbulent  behavior  can  be  considered  fully  established  (AGARD,  1998),  are 
simulated. 


6.1  Test  conditions  and  global  results 

Available  computer  power  made  it  impossible  to  achieve  a  level  of  resolution  correspond¬ 
ing  to  the  one  of  Case  II  (see  Tab.  5.1);  grids  applied  in  present  context  (See  Tab.  6.1)  have 
resolution  roughly  equivalent  to  Case  I  (see  Tab.  5.7). 

It  was  found  that  at  higher  grid  resolution,  results  might  put  in  evidence  the  presence 
of  structures  more  numerous  and  smaller  than  the  ones  to  be  found  by  the  current  grid; 
however,  results  of  §  5.2  compared  to  the  ones  of  §  5.1,  show  that  the  global  information  on 
organized  turbulence  is  well  captured  at  present  grid  resolution;  therefore  it  can  be  assumed 
that  effect  of  Reynolds  number  can  be  reliably  studied  within  this  frame. 


Table  6.1:  Reynolds  effect  study:  Relevant  parameters  for  computational  field  and  grid  sizes. 


Geometry  k  Grid 

El 

E29 

Nx 

Ny 

Nz 

Ax+ 

A  y+ 

A  z+ 

Rer  =  180 

19 

19 

64 

96 

64 

22 

7.85 

1.5 

••  16 

ReT  =  395 

m 

m 

96 

128 

96 

23.27 

6.47 

0.8 

CO 

O 

Rer  =  590 

ta 

EH 

96 

128 

96 

23.25 

6.51 

1.2 

••  44 

It  has  to  be  noticed  that,  even  if  the  computational  field  is  reduced  of  size  with  increasing 
Re,  it  remains  bigger,  by  one  order  of  magnitude,  than  the  longitudinal  macro-scale  in  the 
inertial  layer,  as  yielded  by  DNS  results  (AGARD,  1998);  therefore,  the  corresponding  results 
can  be  considered  to  be  a  representative  sample  of  the  true  flow. 


Detailed  results  of  LES  and  structures  distribution  are  presented  in  Annex  H,  while 
present  discussion  focus  on  the  differences  among  different  Re  s. 

Following  the  conclusions  of  Chapter  5,  Q  criterion  is  applied  to  the  detection  of  struc¬ 
tures;  trigger  level  for  different  cases  is  found  applying  the  criterion  for  decrease  of  PDF  by 
a  decade  (see  §  4.3)  and  the  corresponding  values  are  resumed  in  Tab.  6.2. 

Table  6.2:  Reynolds  effect  study:  Trigger  level  for  different  Rer 


Reynolds  number 

Trigger  level 

Rer  =  180 

72 

Rer  =  395 

185 

Rer  =  590 

271 

Global  values  for  the  two  populations  of  structures  are  compared  in  Tab.  6.3. 

Table  6.3:  Reynolds  effect  study:  Volume,  turbulent  kinetic  energy,  module  of  vorticity  and 
enstrophy  within  structures  versus  corresponding  quantity  within  the  full  computational  field 


Rer 

Volume 

TKE 

V  or.  Mod. 

Ens. 

Num.  Structures 

O 

OO 

t-H 

II 

<u 

ft? 

9.661  % 

14.075  % 

19.893 

% 

25.058  % 

247 

Rer  =  395 

10.926  % 

16.429  % 

25.014 

% 

30.707  % 

323 

Rer  =  590 

11.479  % 

17.413  % 

26.924 

% 

36.854  % 

280 

These  results  point  out  to  a  fundamental  similarity  of  behavior  of  organized  turbulence 
at  all  Reynolds  numbers,  leading  to  the  conclusion  that  what  said  in  Chapter  5  might  be 
considered  valid  for  this  entire  class  of  flows.  It  can  be  remarked  that  contributions  of  orga¬ 
nized  turbulence  in  global  quantities  and  for  the  procedure  here  applied  gradually  increases 
with  Re,  see  at  example  the  field  TKE  in  Tab.  6.3. 

Further  analysis  will  be  performed,  for  the  two  new  Reynolds  number  to  be  considered,  to 
put  in  evidence  the  possible  effects  of  Reynolds  number.  It  must  be  noticed  that,  in  Tab.  6.3, 
the  average  number  of  structures  detected  within  the  computational  field  is  given,  together 
with  information  concerning  global  quantities  (expressed  in  percentage)  of  organized  turbu¬ 
lence;  this  information  is  not  consistent  with  the  others,  since  it  has  not  been  normalized, 
and  it  cannot  be  used  in  order  to  make  comparisons,  for  example  the  domain  size  is  different 
case  by  case. 

However,  this  number  is  given  to  offer  an  indication  of  the  sample  size  to  be  considered. 
It  must  be  kept  in  mind  that  this  number  is  an  averaged  value,  obtained  averaging  data  from, 
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typically,  around  oO  —  60  flow  fields  when  possible  and  in  any  case  no  lower  than  36,  in  order 
to  bring  statistics  to  a  satisfactory  convergence  (see  §  5.1). 

On  the  other  hand,  it  must  be  reminded  that  most  of  the  relevant  features  of  the  pop¬ 
ulation  behavior  (such  as  the  distinguished  peak  in  the  stream-wise  vorticity)  are  already 
discernible  even  over  an  handful  of  gathered  flow-fields  (Giammanco,  2000),  due  to  the  high 
number  of  structures  present  in  one  single  flow  realization. 

Statistics  gathered  on  one  single  snap-shot  are  equivalent  to  an  ensemble  averaging  over 
a  small  population. 
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6.2  Analysis  of  coherent  structures  at  ReT  =  395 

Rer  =  395  corresponds  to  a  bulk  Reynolds  number  of  6875,  high  enough  to  ensure  it 
should  not  anymore  affected  by  the  viscous  effects,  present  in  ReT  =  180.  This  case  has  been 
investigated  in  deep  and  compared  with  reference  DNS  data  to  evaluate  the  quality  of  the 
simulation,  as  shown  in  Annex  H. 

Overall  view  of  vortex  structures  is  presented  in  Fig.  6.1  and  6.2. 


Fig.  6.1:  3-D  overview  of  the  field  of  detected  structures  of  plane  channel  at  Rer  =  395 


Fig.  6.2:  Lateral  overview  of  the  field  of  detected  structures  of  plane  channel  at  ReT  =  395 

Comparison  with  corresponding  views  for  Rer  =  180  (Fig.  5.10  and  5.8)  presents  more 
complex  distribution  of  structures  (Fig.  6.1)  and  increase  of  thickness  of  layer  affected  by  this 
mechanism  (Fig.  6.2).  However  the  overall  features  of  the  field  are  similar. 

The  average  population  is  in  the  order  of  323  structures  (Fig.  6.3)  for  an  average 
value  Volmean  =  4.95e+05  wall  units.  Comparison  with  corresponding  result  at  Rer  =  180 
{Volmean  =  1.48e+05  wall  units)  could  be  considered  an  indication  of  increase  of  typical  sizes 


127 


with  Re.  It  must  remarked  that,  as  discussed  in  §  5.3,  this  value  has  to  be  considered  only 
as  first  crude  indication  of  overall  trend,  while  reliable  estimation  of  size  of  typical  structure 
has  to  be  obtained  by  inertial  moments  of  the  population. 


Fig.  6.3:  Structures  Volume  of  plane  channel  at  ReT  =  395  with  usual  Zone  distinction. 
(Number  of  bins:  8000;  limits  of  histogram:  1.23e+03  •••  6.49e+07) 


Fig.  6.4:  Structures  mass  center  of  plane  channel  at  ReT  =  395  with  usual  zone  distinction. 
(Number  of  bins  200;  limits  of  histogram:  4.7147  •  •  •  790) 

The  distribution  of  structures  mass  center  across  the  height  of  the  channel  (Fig.  6.4) 
shows  a  nearly  uniform  distribution  for  2+  <  200,  indicating  that  the  active  layer  is,  now, 
much  wider  than  one  found  for  previous  case  and  contribution  of  Zone  3  could  no  more  be 
neglected  in  the  analysis  to  be  carried  out  on  the  present  flow:  indeed,  57  %  of  the  structures 
are  now  found  in  Zone  3,  while  Zone  2  now  only  hosts  22%  of  the  structures  found. 
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Advection  velocity  of  mass  center  (Fig.  6.5)  confirms  the  structures  follow  the  average 
advection  of  overall  flow. 


Convection  velocity  along  etream-wise  (wall-unite) 


Fig.  6.5:  Stream-wise  advection  of  structures  center  of  mass  of  plane  channel  at  Rer  =  395. 
(Number  of  bins:  200;  limits  of  histogram:  5.03  •  •  •  20.77) 

Projection  of  structures  on  x,  y,  z  axis  (Fig.  6.6)  shows  sizes  similar  to  ones  found  at 
Rer  =  180  (Fig.  5.15). 

It  can  be  remarked  that  most  probable  value  for  Lx  is  nearly  constant  for  all  the  3  zones 
(Fig.  6.6)  at  Lx  =  100,  while  Ly  and  Lz  increase  with  the  distance  from  the  wall. 

A  more  complete  information  can  be  obtained  from  corresponding  joint  probability  func¬ 
tion  (JPDF)  Lx  -  Ly ;  Ly  -  Lz  and  Lx-Lx.  Distributions  for  the  3  pairing  are  presented  in 
Fig.  6.7  and  Fig.  6.8. 

The  similar  behaviors  of  Zone  2  and  Zone  3  confirms  the  conclusion  that  presence  of 
organized  structures  covers  the  entire  inertial  region  and  corresponding  entities  belong  to  a 
same  phenomenology. 

As  discussed  in  Chapter  5,  the  most  reliable  estimation  of  size  and  orientation  of  the 
structures  is  the  one  obtained  by  the  analysis  of  their  inertial  moments,  which  are  presented 
in  Fig.  6.9  and  Fig.  6.10. 

Present  results  show  patterns  similar  to  ones  discussed  in  §  5.1  and  5.2:  in  Zone  1  most 
of  structures  are  oriented  along  x  direction  (Fig.  6.9  left  column);  in  Zone  2  (Fig.  6.9  right 
column),  turning  of  structures  away  from  x  —  y  plane  increases,  as  well  the  dispersion  of 
structures  in  different  orientations;  this  trend  further  increases  (Fig.  6.10)  in  Zone  3,  making 
it  difficult  to  detect  a  preferred  orientation. 

It  must  be  remarked  that,  since  higher  Reynolds  number  strongly  increases  the  percent- 
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Fig.  6.6:  Structures  statistics  at  ReT  =  395:  length  in  wall  units.  (From  top  left  to 
bottom  right:  number  of  bins:  200,  800,  200;  limits  of  histograms:  46.53  •  •  •  2.23e+03, 
12.92  •  *  •  827.28,  4.22  •  •  •  336.4058) 


age  presence  of  Zone  3  respect  to  Zone  2  and  Zone  1,  being  the  procedure  used  to  produce 
Fig.  6.10  and  Fig.  6.9  the  same  than  the  one  used  to  produce  similar  figures  in  the  previous 
chapters,  the  number  of  structures  present  in  the  Fig.  6.10  appears  extremely  high  compared 
to  Fig.  5.42,  but  this  is  consistent  with  the  relative  presence  of  structures  noted  above. 

Extracting  tilt  and  incidence  angle  for  from  these  data  yields  the  results  presented  in 
Fig.  6.11  to  Fig.  6.14,  where  the  global  normalization  and  the  local  normalization  for  the 
angles  proposed  in  §  5.2  is  applied. 

Fig.  6.11  and  6.13  confirm  that  there  is  no  preferred  tilt  angle;  distribution  of  data  is 
very  similar  for  Zones  1  and  2,  while  variance  of  the  distribution  increases  in  Zone  3,  showing 
that  structures  tend  to  turn  away  from  x  axis  with  increasing  distance  from  wall;  it  can  be 
noted  that  distribution  for  Zone  2  is  similar  to  the  corresponding  one  at  Rer  =  180  but  not 
to  such  an  extent  to  possibly  conclude  of  an  universal  behavior  for  inertial  layer;  for  the  local 
normalization  (conf.  Fig.  5.45  and  Fig.  6.13)  the  PDF  hints  to  a  bigger  dispersion  for  the 
lower  Reynolds. 

Existence  of  an  average  incidence  angle  (Fig.  6.12  and  Fig.  6.13)  shows  the  tendency  of 
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Fig.  6.7:  Structures  statistics:  JPDF  for  structures  extension.  Rer  =  395.  Left  column  Zone 
1,  right  column  Zone  2,  xy- top,  xz-middle  and  yz-bottom. 


structures  to  rotate  away  from  the  x  -  y  plane:  incidence  angle  increases  from  Zone  1  to  3 
to  reach  a  most  probable  state  of  a  “  14  degrees  for  the  intermediate  zone;  again,  it  has  to 
be  remarked  that  the  incidence  angle  for  Zone  2,  here  found,  is  comparable  to  the  one  found 
at  Rer  =  180  both  in  §  5.1  and  §  5.2. 
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Fig.  6.8:  Structures  statistics:  JPDF  for  structures  extension.  Rer  =  395.  Zone  3,  xy-top 
left,  xz-top  right  and  yz-bottom. 


Following  conclusions  reached  in  §  5.1,  the  typical  shape  of  the  structure  is  assumed  to  be 
an  ellipsoid  of  volume  abc,  where  a  >  b  >  c  are  the  semi-length  on  the  3  axis,  whose  values 
are  presented  in  Fig.  6.15  to  Fig.  6.17.  It  is  possible  to  compare  the  volume  of  the  ellipsoid 
with  the  volume  of  the  corresponding  structure,  being  able  finally  to  retrieve  Fig.  6.18. 

According  to  Fig.  6.15  to  6.17,  the  most  probable  state,  for  Zone  2,  would  be  an  ellipsoid 
having  a  =  37.3,  b  =  15.9  and  c  =  11.4  wall  units. 

Values  comparable  to  the  ones  found  in  §  5.1.  It  is  to  be  remarked  that  the  values  2  x  a, 
2x6  are  reasonable  close  to  the  most  probable  state  yielded  by  the  JPDF  Ls  —  Lz  =  75  —  35, 
for  Zone  2  (Fig.  6.7);  this  result  confirms  the  capability  of  present  approach  to  identify  the 
dimensions  of  the  representative  structure. 

Significantly  higher  values  are  found  for  Zone  3:  a  —  57.2,  6  —  25  and  c  =  16.4; 
structures  found  in  Zone  3,  or  in  the  middle  of  the  channel,  are  convected  from  the  walls,  as 
confirmed  by  the  study  of  the  velocity  of  structures  mass  center.  The  fact  that  the  structures 
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Fig.  6.9:  Structures  statistics,  ReT  =  395:  spatial  alignment,  different  pov.  A  Structure 
mean  alignment.  Zone  1  left  column,  Zone  2  right  column 


are  different  in  size  could  hint  to  a  extremely  rapid  growth  process,  but  another  possibility, 
already  addressed  in  §  5.3,  is  that,  close  to  the  wall,  the  structures  are  so  close  to  each  other 
that  the  algorithm  merges  some  of  them  in  some  very  big  structures  with  very  low  probability, 
removing  structures  that  could  possibly  modify  in  a  significantly  way  the  overall  statistics. 
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Fig.  6.10:  Structures  statistics,  Rer  =  395:  spatial  alignment,  different  pov.  A  Structure 
mean  alignment.  Zone  3 


This  fact  alone  should  be  considered  worth  of  further  studies. 

Dynamic  of  organized  vorticity  can  be  assessed  in  Fig.  6.19,  where  specific  vorticity  is 
presented  component  per  component. 

Predominant  clockwise  and  counter-clockwise  rotation  around  x  axis  can  be  noted  in 
the  aforementioned  figure  (top  left),  being  most  evident  for  Zone  3;  this  pattern  of  behavior 
confirms  that  all  structures  present  in  the  inertial  region  belong  to  a  same  family;  different 
behavior  of  Zone  3  with  respect  to  ReT  =  180  case  shows  the  enlarging  of  the  active  region 
and  the  passage  to  fully  developed  turbulent  flow. 

Corresponding  behavior  is  visible  for  wall-normal  specific  vorticity  (bottom  left).  Gaus¬ 
sian  like  behavior  of  specific  vorticity  in  span-wise  direction  (middle  left)  puts  in  evidence 
absence  of  tilt  in  x  -  y  axis,  always  for  average  in  the  whole  channel. 

The  JPDF  of  specific  vorticity  (Fig.  6.20  and  Fig.  6.21)  offer  a  more  complete  under¬ 
standing  of  the  organization  of  vorticity  field;  on  the  basis  of  the  experience  accumulated 


Inerzla  Mom«muums  Driven  Data 


Fig.  6.11:  Structures  statistics  Rer  =  395:  Tilting  angle  extracted  from  inertial  ten¬ 
sor.  (Number  of  bins:  25;  limits  of  histogram:  -90  •  •  •  90) 


Inertia  Mofnertuums  Ortven  Data 


Fig.  6.12:  Structures  statistics  Rer  =  395:  Incidence  angle  extracted  from  inertial  tensor. 
(Number  of  bins:  25;  limits  of  histogram:  —90  •  ■  •  90) 


previous  tests  (see  §  5)  while  for  the  previous  cases,  only  the  behavior  of  the  lower  half  chan¬ 
nel  are  presented  here,  for  clarity  sake,  as  the  joint  PDF  for  the  whole  channel  might  hide 
possible  information. 

The  JPDF  of  specific  vorticity  in  x  and  y  direction  for  all  three  zones  (Fig.  6.20  and 
Fig.  6.21)  show  similar  pattern:  the  overwhelming  majority  of  present  structures  have  a  ro¬ 
tation  along  x;  separation  of  clock-wise  and  counter-clock- wise  rotation  can  be  well  remarked 
in  Zones  2  and  3  and  is  more  clear-cut  than  one  found  for  ReT  =  180.  Values  of  stream-wise 
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Fig.  6.13:,  Structures  statistics  Rer  =  395:  Tilting  angle  extracted  from  inertial  ten¬ 
sor.  (Number  of  bins:25.  Limits  of  histograms  -90  •••  90) 
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Fig.  6.14:  Structures  statistics  Rer  =  395:  Incidence  angle  extracted  from  inertial  ten¬ 
sor.  (Number  of  bins:25.  Limits  of  histograms  —90  •  •  •  90) 


vorticity  are  clustered  around  ±0.10  for  Zone  2  and  ±0.08  for  Zone  3;  values  for  span-wise 
vorticity  are  typically  twice  smaller,  which  is  further  indicating  that  stream-wise  rotation  is 
dominant,  even  if  it  should  not  be  forgotten  that  the  mean  vorticity  has  been  subtracted. 

Similar  behavior  is  evident  for  JPDF  of  specific  vorticity  in  x  and  z  direction  ;  it  has 
to  be  remarked  that,  in  Zone  3,  Fig.  6.21,  wall-normal  specific  vorticity  is  in  the  same  order 
of  magnitude  as  the  one  in  stream-wise  direction,  further  confirming  the  tendency  of  the 
structures  to  turn  away  from  x  —  y  plane. 
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Fig.  6.15:  Structures  statistics  Rer  =  395:  Equivalent  ellipsoid  minor  semi-axis. 


tnerzla  Momentuums  Driven  Data:  Ellipsoid  approach 


Middle  Axl*  (wall-unit*) 


Fig.  6.16:  Structures  statistics  ReT  =  395:  Equivalent  ellipsoid  middle  semi-axis. 


JPDF  of  specific  vorticity  in  y  and  2  direction  shows  the  absence  of  preferred  state  in 
cross-flow  plane. 

6.2.1  Conclusion 

It  can  be  concluded  that  main  findings  for  case  of  ReT  =  180  are  confirmed: 


•  Applying  present  proposal  for  the  choice  of  the  trigger  level  to  case  ReT  =  395  yields 
comparable  values  for  main  physical  parameters  for  organized  turbulence. 


Ifwrzh  Mom#nfuurm  DrK«*o  Data:  EXpaoid  approach 


Fig.  6.17:  Structures  statistics  Rer  =  395:  Equivalent  ellipsoid  major  semi-axis. 


Fig.  6.18:  Structures  statistics  Rer  =  395:  Equivalent  ellipsoid  volume  deviation. 

Most  of  individual  structures  detected  by  the  search  algorithm  have  size  smaller  than 
longitudinal  macro-scale,  even  if  a  minority  of  structures  exists  whose  length  is  compa¬ 
rable  or  bigger  than  this  value. 

Close  to  the  wall  (Zone  1)  structures  are  oriented  along  the  x  axis;  moving  away  from 
the  wall,  turning  away  from  x  —  y  plane  is  observed. 


Results  at  Rer  =  395  are  consistent  with  corresponding  ones  at  ReT  =  180  for  same 
resolution. 


Specific  Vorticity  in  wall-normal  <o2v/uT 
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Fig.  6.19:  Structure  statistics:  Specific  Vorticity  (left)  and  enstrophy (right)  for  x  (top),  y 
(middle)  and  ^(bottom).  Rer  =  395.  (Number  of  bins:  200;  limits  of  histograms  (from  top  to 
bottom),  vorticity:  —0.38  0.32,  —0.46  •••  0.64,  —0.47  •••  0.37;  enstrophy:  0  •••  0.23, 
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Fig.  6.20:  Structures  statistics  lower  half-channel:  JPDF  for  structures  vorticity,  lower  half 
channel.  Left  Zone  1,  right  Zone  2,  ary-top,  xz-middle  and  yz-bottom.  ReT  =  395 
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Fig.  6.21:  Structures  statistics  lower  half-channel:  JPDF  for  structures  vorticity,  lower 
channel.  Zone  3,  xy- top  left,  xz- top  right  and  yz-bottom.Rer  =  395 


6.3  Analysis  of  coherent  structures  at  ReT  =  590 

Complete  results  for  LES  of  eulerian  field  and  organized  structures  are  presented  in 
Annex  I,  while  most  relevant  finding  and  differences  with  results  at  lower  Res  are  reviewed 
below. 

Overall  view  of  vortex  structures  are  presented  in  Fig.  6.22  and  6.23. 
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Fig.  6.22:  3-D  overview  of  the  field  of  detected  structures  of  plane  channel  at  ReT  =  590 
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Fig.  6.23:  Lateral  overview  of  the  field  of  detected  structures  of  plane  channel  at  Rer  =  590 

Comparison  with  corresponding  Fig.  6.1  and  Fig.  6.2  shows  a  remarkable  similarity, 
which  is  a  first  indication  the  region  of  self-similar  behavior  of  turbulence  is  attained;  further 
enlargement  of  the  layer  affected  by  the  process  can  be  remarked. 

The  average  population  is  in  the  order  of  276  structures  (Fig.  6.24)  for  an  average  value 
Volmean  =  9.07e+05  wall  units,  in  line  with  the  tendency  for  volume  of  average  structure 
increases  with  Re,  although  this  value  must  be  considered  as  rough  indication  of  overall 
trend,  as  it  was  already  remarked  discussing  results  for  ReT  =  180  and  ReT  =  395. 
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Fig.  6.24:  Structures  Volume  of  plane  channel  at  ReT  =  590  with  usual  Zone  distinction. 
(Number  of  bins:  2000;  limits  of  histogram:  1.05e+03  •  •  •  l.lle+08) 


Fig.  6.25:  Structures  mass  center  of  plane  channel  at  ReT  =  590  with  usual  Zone  distinction. 
(Number  of  bins:  200;  limits  of  histogram:  3.8  ••  •  1176) 


The  distribution  of  center  of  mass  of  the  structures  across  the  height  of  the  channel 
(Fig.  6.25)  shows  most  of  population  is  now  located  in  region  for  z+  <  300,  confirming  that 
the  active  increases  with  Rer.  The  distribution  appears  nearly  uniform  over  20  <  z+  <  300, 
indicating  all  these  structures  are  part  of  single  phenomenology.  It  has  to  be  remarked  that, 
now  only  18  %  of  structures  belong  to  Zone  2,  while  the  percentage  of  presence  of  structures 
in  Zone  3  is  now  as  high  as  64  %. 

Advection  velocity  of  mass  of  the  structures  (Fig.  6.26)  for  different  zones  is  in  the  same 
order  of  magnitude  as  the  corresponding  average  advection  of  overall  flow,  without  visible 
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Convection  velocity  along  stream-wise  (wall-units) 


Fig.  6.26:  Stream-wise  advection  of  structures  center  of  mass  of  plane  channel  at  ReT  =  590. 
Number  of  bins  200,  limits  3.71  •  •  •  22.11. 


lag. 


Projection  of  structures  on  x,  y  and  z  axis  (Fig.  6.27)  shows  slightly  bigger  sizes  than 
the  ones  found  at  lower  Res. 

It  can  be  remarked  that  most  probable  value  for  Lx  is  nearly  constant  for  all  the  3  zones 
at  Lx  =  100,  while  Ly  and  Lz  increase  with  the  distance  from  the  wall,  to  reach  Ly  =  50  ; 
Lz  —  80  in  Zone  3.  This  trend  further  confirms  the  tendency  of  structures  to  turn  away  from 
x  axis. 

A  more  complete  information  can  be  obtained  from  corresponding  joint  probability  func¬ 
tion  (JPDF)  Lx  —  Ly]  Ly  —  Lz  and  Lx  —  Lx.  Distributions  for  the  3  pairing  are  presented  in 
Fig.  6.28  and  Fig.  6.29. 

Present  results  are  consistent  with  trends  visible  at  lower  Res,  the  most  probable  state 
being:  Lx  —  Lz  =  100  —  40  for  Zone  2  and  Lx  —  Lz  =  100  —  55  for  Zone  3,  which  is  another 
indication  of  turning  away  x  —  y  plane  by  the  structures. 

As  already  stated,  (§  5),  estimation  of  size  and  orientation  of  the  structures  must  be 
obtained  by  the  analysis  of  their  inertial  moments  of  present  population,  which  are  presented 
in  Fig.  6.30  and  6.31. 

Present  results  show  trends  similar  to  ones  observed  for  Rer  =  395  case:  in  Zone  1  most 
of  structures  are  oriented  along  x  direction  (Fig.  6.30  left  column);  in  Zone  2,  turning  of 
structures  away  from  x  —  y  plane  increases  together  with  existence  of  structures  with  different 
orientations  (Fig.  6.30  right  column);  Fig.  6.31  shows  how  most  of  structures  lie  within  Zone 
3;  dispersion  of  data  further  increases,  making  it  difficult  to  recognize  a  preferred  orientation. 
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Extracting  tilt  and  incidence  angle  for  from  these  data  yields  the  results  presented  in 
Fig.  6.32  to' Fig.  6.35. 

8 

Following  the  procedure  introduced  in  §  5.2,  global  normalization  and  local  normalization 
for  the  angles  are  given  as  separate  figures. 

Fig.  6.32  confirms  that  there  is  no  preferred  tilt  angle;  distribution  of  data  is  very  similar 
for  Zones  2  and  3,  variance  of  distribution  being  in  the  same  order  of  magnitude  as  the 
corresponding  one  at  Rer  =  395.  Evolution  of  average  incidence  angle  from  Zone  1  to 
3  shows  the  tendency  of  structures  to  rotate  away  from  the  x  -  y  plane:  incidence  angle 
increases  from  Zone  1  to  3  to  reach  a  most  probable  state  of  a  =  21  degrees  in  Zone  3;  while 
only  7  for  Zone  2.  Again,  the  incidence  angle  for  Zone  2  is  comparable  to  the  one  found  at 
Rer  =  180  and  395,  remembering  conclusions  of  §  5.2. 

The  typical  shape  of  the  structure  is  assumed  to  be  an  ellipsoid  of  volume  §7r abc,  where 
a  >  b  >  c  are  the  semi-length  on  the  3  axis,  whose  values  are  presented  in  Fig.  6.36  to 
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Structures'  Extension  In  span-wise  (wall-units) 
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Fig.  6.28:  Structures  statistics:  JPDF  for  structures  extension.  Rer  =  590.  Left  column 
Zone  1,  right  column  Zone  2,  xy-top,  xx-middle  and  yz-bottom. 


Fig.  6.38.  Corresponding  deviation  is  presented  in  Fig.  6.39. 

Most  probable  values  for  the  three  half-axis,  for  zones  3  are  a  —  54.54  ;  b  —  30.58; 
c  =  17.74.  The  values  are  comparable  to  ones  found  for  Zone  3  at  Rer  =  395,  suggesting  an 
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Fig.  6.29:  Structures  statistics:  JPDF  for  structures  extension.  ReT  =  590.  Zone  3,  xy- top 
left,  xz- top  right  and  yz-bottom. 


universal  phenomenology  for  structures  in  zones  at  high  local  ReT. 

Dynamic  of  organized  vorticity  can  be  assessed  in  Fig.  6.40,  where  specific  vorticity  (see 
§  4)  is  presented  component  per  component. 

Present  trends  are  very  similar  to  ones  remarked  at  at  Rer  =  395:  predominant  clockwise 
and  counter-clockwise  rotation  around  x  axis  can  be  noted  in  the  afore  mentioned  figure(top 
left)  and  is  most  evident  for  Zone  3;  this  pattern  of  behavior  confirms  that  all  structures 
present  in  inertial  region  belong  to  a  same  family.  Corresponding  behavior  is  visible  for 
wall-normal  specific  vorticity  (bottom  left).  Gaussian  like  behavior  of  specific  vorticity  in 
span-wise  direction  (middle  left)  puts  in  evidence  the  lack  of  preferential  orientation  around 
y  axis,  for  whole  channel  averaging. 

JPDF  of  specific  vorticity  in  x  and  y  direction  for  three  zones  (Fig.  6.41  and  Fig.  6.42) 
shows  increasing  organization  moving  away  from  the  wall:  majority  of  present  structures  has 
a  rotation  along  x  ;  separation  of  clock-wise  and  counter-clock-wise  rotation  progresses  from 
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Fig.  6.30:  Structures  statistics,  Rer  =  590:  spatial  alignment,  different  pov.  A  Structure 
mean  alignment.  Zone  1  left  column,  Zone  2  right  column 


Zone  2  to  Zone  3  toward  a  well  defined  preferred  state:  values  of  stream-wise  vorticity  are 
clustered  around  ±0.10  for  Zone  2  and  ±0.08  for  Zone  3;  values  for  span- wise  vorticity  are 
typically  twice  smaller;  all  values  are  very  close  to  the  ones  found  at  ReT  =  395. 


Fig.  6.31:  Structures  statistics,  ReT  =  590:  spatial  alignment,  different  pov.  A  Structure 
mean  alignment.  Zone  3 

Similar  behavior  is  evident  for  JPDF  of  specific  vorticity  in  x  and  2  direction;  it  has  to 
be  remarked  that,  in  Zone  3,  Fig.  6.42,  the  wall-normal  specific  vorticity  is  in  the  same  order 
of  magnitude  as  stream-wise  one,  further  confirming  the  tendency  of  the  structures  to  turn 
away  from  x  —  y  plane. 

JPDF  of  specific  vorticity  in  y  and  z  direction  shows  the  absence  of  preferred  state  in 
cross-flow  plane. 

6.4  Conclusion 

Overall,  present  results  confirm  the  trends  remarked  for  Rer  =  180  and  Rer  =  395, 
and  allow  to  conclude  that,  in  all  three  cases,  organized  turbulence  partakes  the  same  phe¬ 
nomenology. 

It  can  be  concluded  that  the  typical  size  of  structures  and  width  of  layer  affected  by  the 


Inerzla  Momentuums  Driven  Data 


Fig.  6.32:  Structures  statistics  ReT  =  590:  Tilting  angle  extracted  from  inertial  ten- 
sor.(Number  of  bins:25.  Limits  of  histograms  —90  •••  90) 


Inerzta  Momentuuma  Driven  Data 


Fig.  6.33:  Structures  statistics  ReT  =  590:  Incidence  angle  extracted  from  inertial  ten¬ 
sor.  (Number  of  bins:25.  Limits  of  histograms  —90  •••  90) 

phenomenon  seems  to  TO-CHECK  slightly  increase  with  Rer\  pattern  and  values  of  specific 
vorticity  indicated  universal  behavior  of  the  structures  defined  by  present  criterion. 

Results  so  obtained  are  consistent  and  detailed  and,  therefore,  allow  some  general  con¬ 
clusions  which  will  be  presented  in  §  7. 
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Fig.  6.34:  Structures  statistics  Rer  =  590  (Case  I):  Tilting  angle  extracted  from  inertial 
tensor. (Number  of  bins:25.  Limits  of  histograms  -90  •  ■  •  90) 


Inerzia  Momentuums  Driven  Data 


X-Z  Structure  angle  in  degrees 


Fig.  6.35:  Structures  statistics  ReT  =  590:  Incidence  angle  extracted  from  inertial  ten- 
sor.(Number  of  bins:25.  Limits  of  histograms  -90  •  •  •  90) 
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Inwzta  Momwrtuiiim  Drtwi  Data:  ESpaoW  approach 


Fig.  6.36:  Structures  statistics  Rer  =  590:  Equivalent  ellipsoid  minor  semi-axis. 


Fig.  6.37:  Structures  statistics  ReT  =  590:  Equivalent  ellipsoid  middle  semi-axis. 
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Inerzia  Momentuums  Driven  Data:  Ellipsoid  approach 
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Fig.  6.38:  Structures  statistics  Rer  =  590:  Equivalent  ellipsoid  major  semi-axis. 


Inerzia  Momentuums  Driven  Data:  Elfipsokt  approach 
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Fig.  6.39:  Structures  statistics  Rer  -  590:  Equivalent  ellipsoid  volume  deviation. 
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Fig.  6.40:  Structure  statistics:  Specific  Vorticity  (left)  and  enstrophy (right)  for  x  (top),  y 
(middle)  and  z  (bottom).  Rer  =  590.  (Number  of  bins:  200;  limits  of  histograms  (from  top  to 
bottom),  vorticity:  -0.43  •  •  •  0.47,  -0.78  •  •  •  0.62,  -0.49  •  •  •  0.52;  enstrophy:  0  •  •  •  0.3588, 
.29. 
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Fig.  6.42:  Structures  statistics  lower  half-channel:  JPDF  for  structures  vorticity,  lower  half 
channel.  Zone  3,  xy- top  left,  xz- top  right  and  yz-bottom.ReT  =  590 


7.  CONCLUSIONS 


A  search  and  detection  procedure  has  been  developed,  which  allows  to  detect,  define  and 
classify  all  structures  of  organized  vorticity/turbulence  in  3D  flow-fields,  as  well  to  extract  the 
corresponding  relevant  statistics.  The  present  code  is  compatible  with  all  detection  criteria 
proposed  by  the  current  literature  (see  §  4).  Reliability  and  easiness  to  application  lead  to 
the  choice  of  Q  criterion  for  most  present  simulation  and  to  recommend  its  application  for 

future  application. 

It  has  to  be  remarked  that  present  procedure  has  been  developed  and  tested  for  LES 
application,  but  is  compatible  and  is  recommended  for  DNS. 

Concerning  the  study  of  physics  of  coherent  structures,  present  simulations  of  plane 
channel  at  equilibrium  for  3  Reynolds  numbers  allow  to  conclude  that,  in  all  three  cases, 
organized  turbulence  presents  the  same  similar  behavior  and  has  comparable  weight  with 
respect  to  the  overall  turbulent  field.  Therefore,  it  can  be  suggested  that  structures  here 
identified  are  typical  for  attached  wall  layers. 

Application  of  present  detection  and  identification  algorithm  allows  to  conclude  that: 

•  The  procedure  proposed  for  the  choice  of  trigger  level  yields  consistent  values  for  mam 
physical  parameters  for  organized  turbulence  and  can,  therefore,  be  applied  as  general 
rule  to  further  investigations. 

•  Most  of  individual  structures  detected  by  the  search  algorithm  have  sizes  smaller  than 
the  macro-scales  of  the  flow,  even  if  a  minority  of  structures  comparable  or  bigger 
than  this  scale  can  be  found.  It  can  be  remarked  that  volume  of  structure,  expressed 
in  wall  units,  slightly  increases  with  Re,  and  while  the  possibility  that  present  wall- 
unit  adimensionalization  might  not  be  the  most  suited  one  for  present  task  has  been 
addressed,  at  the  moment  it  is  sufficient  to  compare  in  a  meaningful  way  different 

Reynolds. 

•  Close  to  the  wall  (Zone  1)  structures  are  aligned  along  x  axis;  moving  away  from  the 
wall,  a  turning  motion  in  x-y  plane  and  away  from  x-y  plane  is  observed;  typical 
inclination  with  respect  x-y  plane  was  found  to  be  in  order  of  20  degrees. 

It  is  found  that,  refining  the  grid,  global  results  for  the  entire  population  remain  nearly 
constant,  indicating  that  LES  is  adequate  to  recognize  and  assess  the  contribution  of  orga¬ 
nized  structures  to  turbulence;  however,  it  is  also  found  that,  with  increasing  grid  refinement, 
results  indicated  more  and  numerous  and  smaller  structures,  for  similar  volume  occupied  by 
the  entire  population;  therefore  it  is  suggested  that  complete  definition  of  geometry  and  sizes 
of  whole  population  would  be  possible  only  applying  DNS  resolution. 
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APPENDIX 


A.  REt  =  180-CASE  II 

This  appendix  presents  plots  and  information  thought  not  necessary  for  the  main  analysis 
in§  5.2. 


Fig.  A.l:  Anisotropic  part  of  the  resolved  Reynolds  stress  tensor  of  plane  channel  at  Rer  — 
180(Case  II:  refined  grid)  and  total  shear  stress 
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B.  RET  =  180-CASE  IV 


This  appendix  presents  plots  and  information  thought  not  necessary  for  the  main  analysis 
in§  5.4. 


Total  Shear  Stress 


Fig.  B.l:  Anisotropic  part  of  the  resolved  Reynolds  stress  tensor  of  plane  channel  at  ReT  = 
180(Case  IV:  4th  order)  and  total  shear  stress 
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Structures'  Extension  In  span-wise  (wall-units)  Structures’  Extension  In  wall-normal  (wall-units) 

Fig.  B.4:  ReT  =  180(Case  IV:  4th  order),  various  plots 
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'ig.  B.6:  Structures  statistics:  JPDF  for  structures  extension.  Rer  =  180  (Case  IV:  4 
rder) .  Left  column  zone  1,  right  column  zone  2,  xy- top,  xz-middle  and  yz-bottom. 
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Fig.  B.T:  Structure  statistics’:  Specific  Vorticity  (left)  and  enstrophy (right)  for  x  (top),  y 
(middle)  and  z  (bottom). ReT  =  180  (Case  IV:  4t/l  order). 
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D.  REr  =  180-CASE  V-A 


This  appendix  presents  plots  and  information  about  the  FSF  SGS  model  thought  not 
necessary  for  the  main  analysis  in§  5.5. 


Fig.  D.l:  Anisotropic  part  of  the  resolved  Reynolds  stress  tensor  of  plane  channel  at  ReT  = 
180(Case  Y-A:  FSF)  and  total  shear  stress 
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E.  REt  =  1 


This  appen 
not  necessary  fc 


Fig.  E.2:  Mean  stream-wise  velocity  and  rms  of  resolved  turbulence  of  plane  channel  at 
Rer  =  180(Case  V-A:  NO) 


F.  REt  =  180-CASE  V-C 


This  appendix  presents  plots  and  information  about  the  SS  model  thought  not  necessary 
for  the  main  analysis  in§  5.5;  here  the  usual  TKE  budget  is  not  presented  since  the  self  similar 
model  requires  a  different  derivation  that  has  not  been  validated. 


Fig.  F.l:  Mean  stream-wise  velocity  and  rms  of  resolved  turbulence  of  plane  channel  at 
Rer  =  180(Case  V-A:  NO) 
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G.  REt  =  180-CASE  V 


This  appendix  presents  plots  and  informations  are  relative  to  the  comparison  between 
the  simulation  discussed  in  §  5.5,  between  SMA,  FSF,  SS  and  NO-SGS  models. 


Fig.  G.l:  Volume  comparison  between  SGS  models 
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Fig.  G.3:  JPF  for  structure  extensions  between  SGS  models:  Zone  1,  Lx  -  Lv  -II- 
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Fig.  G.5:  JPF  for  structure  extensions  between  SGS  models:  Zone  1 ,  Lx  -  Lz  -II- 
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Structures  Extension  in  wall-normal  (In  wall-units)  Structures  Extension  In  wall-normal  (In  wall-units)  Structures  Extension  In  wall-normal  (In  wall-units) 


Fig.  G.6:  JPF  for  structure  extensions  between  SGS  models:  Zone  1,  Ly  —  Lz  -I- 
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Fig.  G.7:  JPF  for  structure  extensions  between  SGS  models:  Zone  1,  Ly  -  Lz  -II- 
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Fig.  G.8:  JPF  for  structure  extensions  between  SGS  models:  Zone  2,  Lx  —  Ly  -I- 
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Fig.  G.9:  JPF  for  structure  extensions  between  SGS  models:  Zone  2,  Lx  —  Ly  -II 
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Fig.  G.13:  JPF  for  structure  extensions  between  SGS  models:  Zone  2,  Ly  —  Lz  -II- 
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Fig.  G.16:  JPF  for  structure  extensions  between  SGS  models:  Zone  3,  Lx  —  Lz  -I 
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Fig.  G.17:  JPF  for  structure  extensions  between  SGS  models:  Zone  3,  Lx  -  Lz  -II- 
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Fig.  G.18:  JPF  for  structure  extensions  between  SGS  models:  Zone  3,  Ly  —  Lz  -I- 


201 


0  'Cl  20  30  40  50  00  70  00  00  100  “o  10  20  SO  40  00 - 

StnieturM  Extsnalon  In  apan-wlM  (In  waB-unlta)  Structure  Extension  In  apan-wiaa  (In  wall- units) 

Fig.  G.19:  JPF  for  structure  extensions  between  SGS  models:  Zone  3,  Ly  —  Lz  -II- 
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Fig.  G.25:  JPF  for  vorticity  components  between  SGS  models:  Zone  1,  u„  -  u, 
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Fig.  G.29:  JPF  for  vorticity  components  between  SGS  models:  Zone  2,  ux  -  uz  -II- 
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Fig.  G.35:  JPF  for  vorticity  components  between  SGS  models:  Zone  3,  ux  -  uz  -II- 
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Fig.  G.40:  Comparison  for  structures  length  in  wall-normal  between  SGS  models  plus  position 
of  mass  center  in  wall-normal  direction 


Fig.  G.44:  Comparison  for  structures  convection  velocity  in  stream-wise  direction  between 
SGS  models 
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Fig.  G.45:  Comparison  for  structures  convection  velocity  in  wall-normal  direction  between 
SGS  models 
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Fig.  G.50:  Comparison  for  structures  enstrophy  in  span-wise  direction  between  SGS  models 


Fig.  G.51:  Comparison  for  structures  enstrophy  in  wall-normal  direction  between  SGS  models 
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This  appendix  presents  plots  and  information  thought  not  necessary  for  the  main  analysis 
in§  6.2. 
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Fig.  1.3:  Budget  of  resolved  turbulent  kinetic  energy  of  plane  channel  at  Rer  =  590 
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